Synchrony analysis DB-DOS

Libraries

loading dfs

load("data/tidy/tidy_df.RData")
load("data/tidy/tidy_df_pseudo_and_real.RData")
load("data/tidy/tidy_df_pseudo_and_real_per_trial.RData")

Descriptives

Num of dyads

### Number of dyads per group (Pseudo = 0 vs Real = 1)
synchrony_data %>%
  mutate(unique_code = paste0(parent, child)) %>%
  distinct(unique_code, .keep_all = TRUE) %>%
  group_by(category) %>%
  count() %>%
  ungroup()

Missing Channels

### Number of missing channels per group (Pseudo = 0 vs Real = 1 => Category) & Block (1,2,3)
## Values => if 1 = Missing, 0 = not missing
synchrony_data %>%
  mutate(unique_code = paste0(parent, child)) %>%
  select(parent:category) %>%
  pivot_longer(cols = c(s1_d1_hbo:s8_d4_hbo), names_to = "channel", values_to = "values") %>%
  mutate(
    values = case_when(
      is.na(values) ~ "1",
      TRUE ~ "0"
    ),
    values = as.character(values)
  ) %>%
  group_by(category, block, values) %>%
  count() %>%
  ungroup()

Mean Synchrony per Block

### Mean synchrony per block for all Real Dyads
synchrony_data %>%
  filter(category == "Real Dyads") %>%
  group_by(block) %>%
  summarise(mean = round(mean(avg_score, .na.rm = TRUE), 2)) %>%
  ungroup()

Age

### Age
nirs_per_trial_plus_qst %>%
  distinct(record_id, .keep_all = TRUE) %>%
  summarise(
    child_age_mean = round(mean(c_age_y, na.rm = TRUE), 2),
    child_age_sd = round(sd(c_age_y, na.rm = TRUE), 2),
    child_age_min = round(min(c_age_y, na.rm = TRUE), 2),
    child_age_max = round(max(c_age_y, na.rm = TRUE), 2)
  )

Location

### Location
nirs_per_trial_plus_qst %>%
  distinct(record_id, .keep_all = TRUE) %>%
  group_by(site) %>%
  count() %>%
  ungroup()

Sex Parent

### Sex
nirs_per_trial_plus_qst %>%
  distinct(record_id, .keep_all = TRUE) %>%
  group_by(p_sex) %>%
  count() %>%
  ungroup()

Sex Child

### Sex Child
nirs_per_trial_plus_qst %>%
  distinct(record_id, .keep_all = TRUE) %>%
  group_by(c_sex) %>%
  count() %>%
  ungroup()

Race Child

## Race
nirs_per_trial_plus_qst %>%
  distinct(record_id, .keep_all = TRUE) %>%
  group_by(c_race) %>%
  count() %>%
  ungroup() %>%
  filter(c_race != "None") %>%
  mutate(percent = round(n / sum(n) * 100, 2))

Race Parent

## Race Parent

nirs_per_trial_plus_qst %>%
  distinct(record_id, .keep_all = TRUE) %>%
  group_by(p_race) %>%
  count() %>%
  ungroup() %>%
  filter(p_race != "None") %>%
  mutate(percent = round(n / sum(n) * 100, 2))

Latinx Child

## Latinx
nirs_per_trial_plus_qst %>%
  distinct(record_id, .keep_all = TRUE) %>%
  group_by(c_latinx) %>%
  count() %>%
  ungroup() %>%
  filter(c_latinx != "None") %>%
  mutate(percent = round(n / sum(n) * 100, 2))

Latinx Parent

## Latinx Parent
nirs_per_trial_plus_qst %>%
  distinct(record_id, .keep_all = TRUE) %>%
  group_by(p_latinx) %>%
  count() %>%
  ungroup() %>%
  filter(p_latinx != "None") %>%
  mutate(percent = round(n / sum(n) * 100, 2))

Family Income

nirs_per_trial_plus_qst %>%
  distinct(record_id, .keep_all = TRUE) %>%
  group_by(income_house) %>%
  count() %>%
  ungroup() %>%
  filter(income_house != "None") %>%
  mutate(percent = round(n / sum(n) * 100, 2))

Parent Education

nirs_per_trial_plus_qst %>%
  distinct(record_id, .keep_all = TRUE) %>%
  group_by(p_edu) %>%
  count() %>%
  ungroup() %>%
  filter(p_edu != "None") %>%
  mutate(percent = round(n / sum(n) * 100, 2))

Graph of real and Pseudo Dyads

PFC

# synchrony_data$trial <- factor(synchrony_data$trial, levels = c("Trial 1", "Trial 2","Trial 3","Trial 4"))

synchrony_data <- synchrony_data %>%
  mutate(new_block = case_when(
    block == "Block 1" ~ "Free Play Block",
    block == "Block 2" ~ "Stressful Block",
    block == "Block 3" ~ "Recovery Block",
    TRUE ~ "Error"
  ))

synchrony_data$new_block <- factor(synchrony_data$new_block, levels = c("Free Play Block", "Stressful Block", "Recovery Block"))

synchrony_data %>%
  ggplot(., aes(x = category, y = avg_score, fill = category, col = category)) +
  geom_flat_violin(position = position_nudge(x = .2, y = 0), alpha = .2, adjust = 4) +
  geom_point(position = position_jitter(width = .10), size = 1, alpha = 0.4) +
  scale_fill_brewer(palette = "Dark2") +
  scale_colour_brewer(palette = "Dark2") +
  guides(fill = FALSE, col = FALSE) +
  geom_boxplot(aes(x = as.numeric(category) + 0.25, y = avg_score), outlier.shape = NA, alpha = 0.4, width = .1, colour = "BLACK") +
  ylab("Synchrony Values (PFC)") +
  xlab("Dyads Type") +
  stat_summary(fun.data = "mean_cl_normal", geom = "errorbar", width = 1 / 4, color = "black") +
  stat_summary(fun = "mean", geom = "point", color = "black") +
  theme_apa() +
  facet_wrap(~new_block, nrow = 1) +
  theme(text = element_text(size = 30)) +
  scale_x_discrete(labels = function(x) str_wrap(x, width = 8))

## Warning: The `<scale>` argument of `guides()` cannot be `FALSE`. Use "none" instead as
## of ggplot2 3.3.4.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

## Warning: Using the `size` aesthetic with geom_polygon was deprecated in ggplot2 3.4.0.
## ℹ Please use the `linewidth` aesthetic instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

ggsave("figures/aim_0_real_pseudo_dyads_avg-PFC.png", width = 14, height = 10)

DLPFC

synchrony_data %>%
  ggplot(., aes(x = category, y = DLPFC_avg_score, fill = category, col = category)) +
  geom_flat_violin(position = position_nudge(x = .2, y = 0), alpha = .2, adjust = 4) +
  geom_point(position = position_jitter(width = .10), size = 1, alpha = 0.4) +
  scale_fill_brewer(palette = "Dark2") +
  scale_colour_brewer(palette = "Dark2") +
  guides(fill = FALSE, col = FALSE) +
  geom_boxplot(aes(x = as.numeric(category) + 0.25, y = DLPFC_avg_score), outlier.shape = NA, alpha = 0.4, width = .1, colour = "BLACK") +
  ylab("Synchrony Values (DLPFC)") +
  xlab("Dyads Type") +
  stat_summary(fun.data = "mean_cl_normal", geom = "errorbar", width = 1 / 4, color = "black") +
  stat_summary(fun = "mean", geom = "point", color = "black") +
  theme_apa() +
  facet_wrap(~new_block, nrow = 1) +
  theme(text = element_text(size = 20)) +
  scale_x_discrete(labels = function(x) str_wrap(x, width = 8))

## Warning: Removed 14 rows containing non-finite outside the scale range
## (`stat_ydensity()`).

## Warning: Removed 14 rows containing non-finite outside the scale range
## (`stat_boxplot()`).

## Warning: Removed 14 rows containing non-finite outside the scale range
## (`stat_summary()`).
## Removed 14 rows containing non-finite outside the scale range
## (`stat_summary()`).

## Warning: Removed 14 rows containing missing values or values outside the scale range
## (`geom_point()`).

ggsave("figures/aim_0_real_pseudo_dyads_avg-DLPFC.png", width = 14, height = 10)

## Warning: Removed 14 rows containing non-finite outside the scale range
## (`stat_ydensity()`).

## Warning: Removed 14 rows containing non-finite outside the scale range
## (`stat_boxplot()`).

## Warning: Removed 14 rows containing non-finite outside the scale range
## (`stat_summary()`).
## Removed 14 rows containing non-finite outside the scale range
## (`stat_summary()`).

## Warning: Removed 14 rows containing missing values or values outside the scale range
## (`geom_point()`).

VLPFC

synchrony_data %>%
  ggplot(., aes(x = category, y = VLPFC_avg_score, fill = category, col = category)) +
  geom_flat_violin(position = position_nudge(x = .2, y = 0), alpha = .2, adjust = 4) +
  geom_point(position = position_jitter(width = .10), size = 1, alpha = 0.4) +
  scale_fill_brewer(palette = "Dark2") +
  scale_colour_brewer(palette = "Dark2") +
  guides(fill = FALSE, col = FALSE) +
  geom_boxplot(aes(x = as.numeric(category) + 0.25, y = VLPFC_avg_score), outlier.shape = NA, alpha = 0.4, width = .1, colour = "BLACK") +
  ylab("Synchrony Values (VLPFC)") +
  xlab("Dyads Type") +
  stat_summary(fun.data = "mean_cl_normal", geom = "errorbar", width = 1 / 4, color = "black") +
  stat_summary(fun = "mean", geom = "point", color = "black") +
  theme_apa() +
  facet_wrap(~new_block, nrow = 1) +
  theme(text = element_text(size = 20)) +
  scale_x_discrete(labels = function(x) str_wrap(x, width = 8))

ggsave("figures/aim_0_real_pseudo_dyads_avg-VLPFC.png", width = 14, height = 10)

Model Testing Pseudo and real dyads

nirs_per_trial_plus_qst$block <- factor(nirs_per_trial_plus_qst$block, levels = c("Block 2", "Block 1", "Block 3"))

nirs_per_trial_plus_qst <- nirs_per_trial_plus_qst %>%
  mutate(new_id = paste0(parent, child))


model1 <- lmer(
  avg_score ~ category +
    (1 + block | record_id) + (1 | trial),
  data = synchrony_data_per_trial,
  control = lmerControl(
    optimizer = "bobyqa",
    optCtrl = list(maxfun = 5e5)
  )
)

summary(model1)

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: avg_score ~ category + (1 + block | record_id) + (1 | trial)
##    Data: synchrony_data_per_trial
## Control: lmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 5e+05))
## 
## REML criterion at convergence: -557606.4
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -5.5789 -0.6662 -0.0163  0.6408  6.3306 
## 
## Random effects:
##  Groups    Name         Variance  Std.Dev.  Corr       
##  record_id (Intercept)  5.161e-05 0.0071843            
##            blockBlock 2 6.553e-04 0.0255983 -0.87      
##            blockBlock 3 3.039e-05 0.0055124 -0.35  0.09
##  trial     (Intercept)  9.209e-07 0.0009596            
##  Residual               1.128e-03 0.0335837            
## Number of obs: 141451, groups:  record_id, 109; trial, 4
## 
## Fixed effects:
##                        Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)           3.473e-01  1.084e-03  7.271e+01 320.401  < 2e-16 ***
## categoryPseudo Dyads -4.050e-03  9.347e-04  1.411e+05  -4.333 1.47e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## ctgryPsdDyd -0.854

report(model1)

## We fitted a linear mixed model (estimated using REML and BOBYQA optimizer) to
## predict avg_score with category (formula: avg_score ~ category). The model
## included block as random effects (formula: list(~1 + block | record_id, ~1 |
## trial)). The model's total explanatory power is weak (conditional R2 = 0.04)
## and the part related to the fixed effects alone (marginal R2) is of 1.27e-04.
## The model's intercept, corresponding to category = Real Dyads, is at 0.35 (95%
## CI [0.35, 0.35], t(141441) = 320.40, p < .001). Within this model:
## 
##   - The effect of category [Pseudo Dyads] is statistically significant and
## negative (beta = -4.05e-03, 95% CI [-5.88e-03, -2.22e-03], t(141441) = -4.33, p
## < .001; Std. beta = -0.11, 95% CI [-0.16, -0.06])
## 
## Standardized parameters were obtained by fitting the model on a standardized
## version of the dataset. 95% Confidence Intervals (CIs) and p-values were
## computed using a Wald t-distribution approximation.

Table Summary

modelsummary_1 <- list(
  "Full Model" = model1
)

msummary(modelsummary_1,
  output = "tables/aim_0.docx",
  stars = TRUE, align = "lc", metrics = c("RMSE", "R2", "AIC", "BIC", "Log.Lik.", "F")
)

Sum of blank channels

synchrony_data %>%
  filter(category == "Real Dyads") %>%
  dplyr::mutate(across(c("s1_d1_hbo":"s8_d4_hbo"), as.numeric)) %>%
  rowwise() %>%
  mutate(num_na = sum(is.na(c_across(s1_d1_hbo:s8_d4_hbo)))) %>%
  group_by(block) %>%
  summarise(num_na = sum(num_na, na.rm = TRUE)) %>%
  ungroup()

Adding Lab ID

1 = WUSTL 0 = PSU

# nirs_per_trial_plus_qst$block <- factor(nirs_per_trial_plus_qst$block, levels = #c("Block 2", "Block 1", "Block 3"))
nirs_per_trial_plus_qst <- nirs_per_trial_plus_qst %>%
  mutate(
    block_2 = case_when(
      block == "Block 1" ~ 1,
      block == "Block 2" ~ 2,
      block == "Block 3" ~ 3
    ),
    block_sq = block_2^2
  )


nirs_per_trial_plus_qst <- nirs_per_trial_plus_qst %>%
  mutate(lab = case_when(
    record_id >= 1200 ~ 0,
    TRUE ~ 1
  ))

Breaking apart datasets

block1 <- nirs_per_trial_plus_qst %>%
  filter(block == "Block 1")

block2 <- nirs_per_trial_plus_qst %>%
  filter(block == "Block 2")

block3 <- nirs_per_trial_plus_qst %>%
  filter(block == "Block 3")

Correlations

p-values(0, 0.001, 0.01, 0.05, 0.1, 1) <=> symbols(“”, “”, “”, “.”, ” “)

correlation_df <- block1 %>%
  select(stai_trait, stai_state, cbq_effortful_control, cbq_fear, anxiety_c)

chart.Correlation(correlation_df, histogram = TRUE, method = "pearson")

Correlations-block 1 neuro

correlation_df <- block1 %>%
  select(stai_trait, stai_state, cbq_effortful_control, cbq_fear, anxiety_c, avg_score, DLPFC_avg_score, VLPFC_avg_score)

chart.Correlation(correlation_df, histogram = TRUE, method = "pearson")

Correlations-block 2 neuro

correlation_df <- block2 %>%
  select(stai_trait, stai_state, cbq_effortful_control, cbq_fear, anxiety_c, avg_score, DLPFC_avg_score, VLPFC_avg_score)

chart.Correlation(correlation_df, histogram = TRUE, method = "pearson")

Correlations-block 3 neuro

correlation_df <- block3 %>%
  select(stai_trait, stai_state, cbq_effortful_control, cbq_fear, anxiety_c, avg_score, DLPFC_avg_score, VLPFC_avg_score)

chart.Correlation(correlation_df, histogram = TRUE, method = "pearson")

Multilevel Correlations

mutilevel_cor_data <- nirs_per_trial_plus_qst %>%
  select(block, trial, cbq_fear, cbq_effortful_control, anxiety_c, stai_state, stai_trait, avg_score, DLPFC_avg_score, VLPFC_avg_score) %>%
  mutate(
    block = as_factor(block),
    trial = as_factor(trial)
  )

correlation::correlation(mutilevel_cor_data, multilevel = TRUE)

## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')

Lab Comparisons

there is no significant interaction between lab and blocks.

model0 <- lmer(avg_score ~ site + (1 + block | record_id) + (1 | trial),
  data = nirs_per_trial_plus_qst,
  control = lmerControl(
    optimizer = "bobyqa",
    optCtrl = list(maxfun = 5e5)
  )
)

## boundary (singular) fit: see help('isSingular')

summary(model0)

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: avg_score ~ site + (1 + block | record_id) + (1 | trial)
##    Data: nirs_per_trial_plus_qst
## Control: lmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 5e+05))
## 
## REML criterion at convergence: -5024.1
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.9120 -0.6408 -0.0343  0.6387  3.5459 
## 
## Random effects:
##  Groups    Name         Variance  Std.Dev. Corr       
##  record_id (Intercept)  0.0004310 0.02076             
##            blockBlock 1 0.0007491 0.02737  -1.00      
##            blockBlock 3 0.0005875 0.02424  -0.98  0.99
##  trial     (Intercept)  0.0000000 0.00000             
##  Residual               0.0011146 0.03339             
## Number of obs: 1303, groups:  record_id, 109; trial, 4
## 
## Fixed effects:
##               Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)   0.348847   0.001449 200.051826 240.747   <2e-16 ***
## siteWUSTL    -0.001632   0.001950 198.725133  -0.837    0.404    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##           (Intr)
## siteWUSTL -0.743
## optimizer (bobyqa) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')

report(model0)

## boundary (singular) fit: see help('isSingular')

## Random effect variances not available. Returned R2 does not account for random effects.

## boundary (singular) fit: see help('isSingular')

## Random effect variances not available. Returned R2 does not account for random effects.

## We fitted a linear mixed model (estimated using REML and BOBYQA optimizer) to
## predict avg_score with site (formula: avg_score ~ site). The model included
## block as random effects (formula: list(~1 + block | record_id, ~1 | trial)).
## The model's explanatory power related to the fixed effects alone (marginal R2)
## is 5.91e-04. The model's intercept, corresponding to site = PSU, is at 0.35
## (95% CI [0.35, 0.35], t(1293) = 240.75, p < .001). Within this model:
## 
##   - The effect of site [WUSTL] is statistically non-significant and negative
## (beta = -1.63e-03, 95% CI [-5.46e-03, 2.19e-03], t(1293) = -0.84, p = 0.403;
## Std. beta = -0.05, 95% CI [-0.15, 0.06])
## 
## Standardized parameters were obtained by fitting the model on a standardized
## version of the dataset. 95% Confidence Intervals (CIs) and p-values were
## computed using a Wald t-distribution approximation.

model00 <- lmer(DLPFC_avg_score ~ site + (1 + block | record_id) + (1 | trial),
  data = nirs_per_trial_plus_qst,
  control = lmerControl(
    optimizer = "bobyqa",
    optCtrl = list(maxfun = 5e5)
  )
)

## boundary (singular) fit: see help('isSingular')

summary(model00)

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: DLPFC_avg_score ~ site + (1 + block | record_id) + (1 | trial)
##    Data: nirs_per_trial_plus_qst
## Control: lmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 5e+05))
## 
## REML criterion at convergence: -3706.7
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.2474 -0.6846  0.0063  0.6553  3.5794 
## 
## Random effects:
##  Groups    Name         Variance  Std.Dev. Corr       
##  record_id (Intercept)  0.0006577 0.02565             
##            blockBlock 1 0.0007975 0.02824  -1.00      
##            blockBlock 3 0.0012254 0.03501  -1.00  1.00
##  trial     (Intercept)  0.0000000 0.00000             
##  Residual               0.0028238 0.05314             
## Number of obs: 1254, groups:  record_id, 109; trial, 4
## 
## Fixed effects:
##               Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)   0.350059   0.002309 611.688877 151.607   <2e-16 ***
## siteWUSTL    -0.001663   0.003090 606.170555  -0.538    0.591    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##           (Intr)
## siteWUSTL -0.747
## optimizer (bobyqa) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')

report(model00)

## boundary (singular) fit: see help('isSingular')

## Random effect variances not available. Returned R2 does not account for random effects.

## boundary (singular) fit: see help('isSingular')

## Random effect variances not available. Returned R2 does not account for random effects.

## We fitted a linear mixed model (estimated using REML and BOBYQA optimizer) to
## predict DLPFC_avg_score with site (formula: DLPFC_avg_score ~ site). The model
## included block as random effects (formula: list(~1 + block | record_id, ~1 |
## trial)). The model's explanatory power related to the fixed effects alone
## (marginal R2) is 2.42e-04. The model's intercept, corresponding to site = PSU,
## is at 0.35 (95% CI [0.35, 0.35], t(1244) = 151.61, p < .001). Within this
## model:
## 
##   - The effect of site [WUSTL] is statistically non-significant and negative
## (beta = -1.66e-03, 95% CI [-7.72e-03, 4.40e-03], t(1244) = -0.54, p = 0.591;
## Std. beta = -0.03, 95% CI [-0.14, 0.08])
## 
## Standardized parameters were obtained by fitting the model on a standardized
## version of the dataset. 95% Confidence Intervals (CIs) and p-values were
## computed using a Wald t-distribution approximation.

model000 <- lmer(avg_score ~ site + (1 + block | record_id) + (1 | trial),
  data = nirs_per_trial_plus_qst,
  control = lmerControl(
    optimizer = "bobyqa",
    optCtrl = list(maxfun = 5e5)
  )
)

## boundary (singular) fit: see help('isSingular')

summary(model000)

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: avg_score ~ site + (1 + block | record_id) + (1 | trial)
##    Data: nirs_per_trial_plus_qst
## Control: lmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 5e+05))
## 
## REML criterion at convergence: -5024.1
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.9120 -0.6408 -0.0343  0.6387  3.5459 
## 
## Random effects:
##  Groups    Name         Variance  Std.Dev. Corr       
##  record_id (Intercept)  0.0004310 0.02076             
##            blockBlock 1 0.0007491 0.02737  -1.00      
##            blockBlock 3 0.0005875 0.02424  -0.98  0.99
##  trial     (Intercept)  0.0000000 0.00000             
##  Residual               0.0011146 0.03339             
## Number of obs: 1303, groups:  record_id, 109; trial, 4
## 
## Fixed effects:
##               Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)   0.348847   0.001449 200.051826 240.747   <2e-16 ***
## siteWUSTL    -0.001632   0.001950 198.725133  -0.837    0.404    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##           (Intr)
## siteWUSTL -0.743
## optimizer (bobyqa) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')

report(model000)

## boundary (singular) fit: see help('isSingular')

## Random effect variances not available. Returned R2 does not account for random effects.

## boundary (singular) fit: see help('isSingular')

## Random effect variances not available. Returned R2 does not account for random effects.

## We fitted a linear mixed model (estimated using REML and BOBYQA optimizer) to
## predict avg_score with site (formula: avg_score ~ site). The model included
## block as random effects (formula: list(~1 + block | record_id, ~1 | trial)).
## The model's explanatory power related to the fixed effects alone (marginal R2)
## is 5.91e-04. The model's intercept, corresponding to site = PSU, is at 0.35
## (95% CI [0.35, 0.35], t(1293) = 240.75, p < .001). Within this model:
## 
##   - The effect of site [WUSTL] is statistically non-significant and negative
## (beta = -1.63e-03, 95% CI [-5.46e-03, 2.19e-03], t(1293) = -0.84, p = 0.403;
## Std. beta = -0.05, 95% CI [-0.15, 0.06])
## 
## Standardized parameters were obtained by fitting the model on a standardized
## version of the dataset. 95% Confidence Intervals (CIs) and p-values were
## computed using a Wald t-distribution approximation.

1. Aim 1: To examine parent-child neural synchrony using the DBDOS BioSync version

Main Models

PFC

nirs_per_trial_plus_qst$block <- factor(nirs_per_trial_plus_qst$block, levels = c("Block 2", "Block 1", "Block 3"))

model1 <- lmer(
  avg_score ~
    c_age_y_centered + block + c_sex + (1 + block | record_id) + (1 | trial),
  data = nirs_per_trial_plus_qst,
  control = lmerControl(
    optimizer = "bobyqa",
    optCtrl = list(maxfun = 5e5)
  )
)

## boundary (singular) fit: see help('isSingular')

summary(model1)

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: 
## avg_score ~ c_age_y_centered + block + c_sex + (1 + block | record_id) +  
##     (1 | trial)
##    Data: nirs_per_trial_plus_qst
## Control: lmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 5e+05))
## 
## REML criterion at convergence: -5052
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.8782 -0.6598 -0.0480  0.6164  3.7013 
## 
## Random effects:
##  Groups    Name         Variance  Std.Dev. Corr       
##  record_id (Intercept)  0.0001765 0.01329             
##            blockBlock 1 0.0002932 0.01712  -1.00      
##            blockBlock 3 0.0001386 0.01177  -1.00  1.00
##  trial     (Intercept)  0.0000000 0.00000             
##  Residual               0.0011160 0.03341             
## Number of obs: 1303, groups:  record_id, 109; trial, 4
## 
## Fixed effects:
##                    Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)       3.638e-01  2.276e-03  1.539e+02 159.896  < 2e-16 ***
## c_age_y_centered  5.723e-04  8.258e-04  3.631e+02   0.693    0.489    
## blockBlock 1     -2.129e-02  2.797e-03  1.178e+02  -7.614 7.26e-12 ***
## blockBlock 3     -2.151e-02  2.531e-03  1.498e+02  -8.497 1.80e-14 ***
## c_sexF            1.675e-04  1.981e-03  3.640e+02   0.085    0.933    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) c_g_y_ blckB1 blckB3
## c_g_y_cntrd  0.093                     
## blockBlock1 -0.731  0.000              
## blockBlock3 -0.692 -0.001  0.623       
## c_sexF      -0.439 -0.213  0.002 -0.002
## optimizer (bobyqa) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')

report(model1)

## boundary (singular) fit: see help('isSingular')

## Random effect variances not available. Returned R2 does not account for random effects.

## boundary (singular) fit: see help('isSingular')

## Random effect variances not available. Returned R2 does not account for random effects.

## We fitted a linear mixed model (estimated using REML and BOBYQA optimizer) to
## predict avg_score with c_age_y_centered, block and c_sex (formula: avg_score ~
## c_age_y_centered + block + c_sex). The model included block as random effects
## (formula: list(~1 + block | record_id, ~1 | trial)). The model's explanatory
## power related to the fixed effects alone (marginal R2) is 0.08. The model's
## intercept, corresponding to c_age_y_centered = 0, block = Block 2 and c_sex =
## M, is at 0.36 (95% CI [0.36, 0.37], t(1290) = 159.90, p < .001). Within this
## model:
## 
##   - The effect of c age y centered is statistically non-significant and positive
## (beta = 5.72e-04, 95% CI [-1.05e-03, 2.19e-03], t(1290) = 0.69, p = 0.488; Std.
## beta = 0.02, 95% CI [-0.04, 0.07])
##   - The effect of block [Block 1] is statistically significant and negative (beta
## = -0.02, 95% CI [-0.03, -0.02], t(1290) = -7.61, p < .001; Std. beta = -0.60,
## 95% CI [-0.75, -0.44])
##   - The effect of block [Block 3] is statistically significant and negative (beta
## = -0.02, 95% CI [-0.03, -0.02], t(1290) = -8.50, p < .001; Std. beta = -0.60,
## 95% CI [-0.74, -0.46])
##   - The effect of c sex [F] is statistically non-significant and positive (beta =
## 1.68e-04, 95% CI [-3.72e-03, 4.05e-03], t(1290) = 0.08, p = 0.933; Std. beta =
## 4.68e-03, 95% CI [-0.10, 0.11])
## 
## Standardized parameters were obtained by fitting the model on a standardized
## version of the dataset. 95% Confidence Intervals (CIs) and p-values were
## computed using a Wald t-distribution approximation.

performance::check_outliers(model1)

## OK: No outliers detected.
## - Based on the following method and threshold: cook (0.8).
## - For variable: (Whole model)

DLPFC

nirs_per_trial_plus_qst$block <- factor(nirs_per_trial_plus_qst$block, levels = c("Block 2", "Block 1", "Block 3"))

model2 <- lmer(
  DLPFC_avg_score ~
    c_age_y_centered + block + c_sex + (1 + block | record_id) + (1 | trial),
  data = nirs_per_trial_plus_qst,
  control = lmerControl(
    optimizer = "bobyqa",
    optCtrl = list(maxfun = 5e5)
  )
)

## boundary (singular) fit: see help('isSingular')

summary(model2)

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: DLPFC_avg_score ~ c_age_y_centered + block + c_sex + (1 + block |  
##     record_id) + (1 | trial)
##    Data: nirs_per_trial_plus_qst
## Control: lmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 5e+05))
## 
## REML criterion at convergence: -3711.4
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.2899 -0.6839 -0.0244  0.6429  3.6518 
## 
## Random effects:
##  Groups    Name         Variance  Std.Dev. Corr       
##  record_id (Intercept)  0.0003854 0.01963             
##            blockBlock 1 0.0003919 0.01980  -1.00      
##            blockBlock 3 0.0006443 0.02538  -1.00  1.00
##  trial     (Intercept)  0.0000000 0.00000             
##  Residual               0.0028141 0.05305             
## Number of obs: 1254, groups:  record_id, 109; trial, 4
## 
## Fixed effects:
##                    Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)       3.662e-01  3.595e-03  1.488e+02 101.867  < 2e-16 ***
## c_age_y_centered  2.704e-03  1.321e-03  4.749e+02   2.047   0.0412 *  
## blockBlock 1     -2.102e-02  4.145e-03  1.387e+02  -5.071 1.25e-06 ***
## blockBlock 3     -2.447e-02  4.405e-03  1.111e+02  -5.555 1.93e-07 ***
## c_sexF           -6.878e-04  3.177e-03  4.744e+02  -0.217   0.8287    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) c_g_y_ blckB1 blckB3
## c_g_y_cntrd  0.107                     
## blockBlock1 -0.699 -0.002              
## blockBlock3 -0.718 -0.001  0.626       
## c_sexF      -0.446 -0.242  0.004 -0.003
## optimizer (bobyqa) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')

report(model2)

## boundary (singular) fit: see help('isSingular')

## Random effect variances not available. Returned R2 does not account for random effects.

## boundary (singular) fit: see help('isSingular')

## Random effect variances not available. Returned R2 does not account for random effects.

## We fitted a linear mixed model (estimated using REML and BOBYQA optimizer) to
## predict DLPFC_avg_score with c_age_y_centered, block and c_sex (formula:
## DLPFC_avg_score ~ c_age_y_centered + block + c_sex). The model included block
## as random effects (formula: list(~1 + block | record_id, ~1 | trial)). The
## model's explanatory power related to the fixed effects alone (marginal R2) is
## 0.04. The model's intercept, corresponding to c_age_y_centered = 0, block =
## Block 2 and c_sex = M, is at 0.37 (95% CI [0.36, 0.37], t(1241) = 101.87, p <
## .001). Within this model:
## 
##   - The effect of c age y centered is statistically significant and positive
## (beta = 2.70e-03, 95% CI [1.12e-04, 5.30e-03], t(1241) = 2.05, p = 0.041; Std.
## beta = 0.06, 95% CI [2.43e-03, 0.11])
##   - The effect of block [Block 1] is statistically significant and negative (beta
## = -0.02, 95% CI [-0.03, -0.01], t(1241) = -5.07, p < .001; Std. beta = -0.38,
## 95% CI [-0.53, -0.23])
##   - The effect of block [Block 3] is statistically significant and negative (beta
## = -0.02, 95% CI [-0.03, -0.02], t(1241) = -5.55, p < .001; Std. beta = -0.44,
## 95% CI [-0.60, -0.29])
##   - The effect of c sex [F] is statistically non-significant and negative (beta =
## -6.88e-04, 95% CI [-6.92e-03, 5.54e-03], t(1241) = -0.22, p = 0.829; Std. beta
## = -0.01, 95% CI [-0.12, 0.10])
## 
## Standardized parameters were obtained by fitting the model on a standardized
## version of the dataset. 95% Confidence Intervals (CIs) and p-values were
## computed using a Wald t-distribution approximation.

performance::check_outliers(model2)

## OK: No outliers detected.
## - Based on the following method and threshold: cook (0.8).
## - For variable: (Whole model)

VLPFC

nirs_per_trial_plus_qst$block <- factor(nirs_per_trial_plus_qst$block, levels = c("Block 2", "Block 1", "Block 3"))

model3 <- lmer(
  VLPFC_avg_score ~
    c_age_y_centered + block + c_sex + (1 + block | record_id) + (1 | trial),
  data = nirs_per_trial_plus_qst,
  control = lmerControl(
    optimizer = "bobyqa",
    optCtrl = list(maxfun = 5e5)
  )
)

## boundary (singular) fit: see help('isSingular')

summary(model3)

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: VLPFC_avg_score ~ c_age_y_centered + block + c_sex + (1 + block |  
##     record_id) + (1 | trial)
##    Data: nirs_per_trial_plus_qst
## Control: lmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 5e+05))
## 
## REML criterion at convergence: -4725.7
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.7868 -0.6707  0.0030  0.6457  4.2572 
## 
## Random effects:
##  Groups    Name         Variance  Std.Dev. Corr       
##  record_id (Intercept)  9.633e-05 0.009815            
##            blockBlock 1 1.230e-04 0.011090 -0.99      
##            blockBlock 3 5.317e-05 0.007292 -0.93  0.97
##  trial     (Intercept)  0.000e+00 0.000000            
##  Residual               1.452e-03 0.038106            
## Number of obs: 1300, groups:  record_id, 109; trial, 4
## 
## Fixed effects:
##                    Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)       3.626e-01  2.360e-03  1.531e+02 153.688  < 2e-16 ***
## c_age_y_centered -5.506e-04  9.571e-04  1.280e+02  -0.575    0.566    
## blockBlock 1     -2.173e-02  2.798e-03  1.279e+02  -7.767 2.28e-12 ***
## blockBlock 3     -2.094e-02  2.683e-03  1.666e+02  -7.805 6.24e-13 ***
## c_sexF           -1.290e-05  2.295e-03  1.284e+02  -0.006    0.996    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) c_g_y_ blckB1 blckB3
## c_g_y_cntrd  0.105                     
## blockBlock1 -0.658 -0.002              
## blockBlock3 -0.624 -0.002  0.543       
## c_sexF      -0.488 -0.212  0.000 -0.005
## optimizer (bobyqa) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')

report(model3)

## boundary (singular) fit: see help('isSingular')

## Random effect variances not available. Returned R2 does not account for random effects.

## boundary (singular) fit: see help('isSingular')

## Random effect variances not available. Returned R2 does not account for random effects.

## We fitted a linear mixed model (estimated using REML and BOBYQA optimizer) to
## predict VLPFC_avg_score with c_age_y_centered, block and c_sex (formula:
## VLPFC_avg_score ~ c_age_y_centered + block + c_sex). The model included block
## as random effects (formula: list(~1 + block | record_id, ~1 | trial)). The
## model's explanatory power related to the fixed effects alone (marginal R2) is
## 0.07. The model's intercept, corresponding to c_age_y_centered = 0, block =
## Block 2 and c_sex = M, is at 0.36 (95% CI [0.36, 0.37], t(1287) = 153.69, p <
## .001). Within this model:
## 
##   - The effect of c age y centered is statistically non-significant and negative
## (beta = -5.51e-04, 95% CI [-2.43e-03, 1.33e-03], t(1287) = -0.58, p = 0.565;
## Std. beta = -0.02, 95% CI [-0.07, 0.04])
##   - The effect of block [Block 1] is statistically significant and negative (beta
## = -0.02, 95% CI [-0.03, -0.02], t(1287) = -7.77, p < .001; Std. beta = -0.55,
## 95% CI [-0.68, -0.41])
##   - The effect of block [Block 3] is statistically significant and negative (beta
## = -0.02, 95% CI [-0.03, -0.02], t(1287) = -7.81, p < .001; Std. beta = -0.53,
## 95% CI [-0.66, -0.39])
##   - The effect of c sex [F] is statistically non-significant and negative (beta =
## -1.29e-05, 95% CI [-4.52e-03, 4.49e-03], t(1287) = -5.62e-03, p = 0.996; Std.
## beta = -3.24e-04, 95% CI [-0.11, 0.11])
## 
## Standardized parameters were obtained by fitting the model on a standardized
## version of the dataset. 95% Confidence Intervals (CIs) and p-values were
## computed using a Wald t-distribution approximation.

performance::check_outliers(model3)

## OK: No outliers detected.
## - Based on the following method and threshold: cook (0.8).
## - For variable: (Whole model)

PFC Graph

synchrony_data$block <- factor(synchrony_data$block, levels = c("Block 1", "Block 2", "Block 3"))

synchrony_data %>%
  filter(category == "Real Dyads") %>%
  ggplot(., aes(x = new_block, y = avg_score, fill = block, col = block)) +
  scale_color_manual(
    name = "Block Type",
    values = c("#C85729", "#127088", "#92874B")
  ) +
  scale_fill_manual(
    name = "Block Type",
    values = c("#C85729", "#127088", "#92874B")
  ) +
  guides(fill = FALSE, col = FALSE) +
  geom_rain(alpha = .5, rain.side = "f2x2", id.long.var = "record_id") +
  ylab("Neural Synchrony Values (PFC)") +
  xlab("Block Type") +
  stat_summary(fun.data = "mean_cl_normal", geom = "errorbar", width = 1 / 4, color = "black") +
  stat_summary(fun = "mean", geom = "point", color = "black") +
  theme_apa() +
  theme(text = element_text(size = 15))

## Warning: Duplicated aesthetics after name standardisation: alpha

## Warning in setup_data(...): geom_paired_raincloud is only useful for visualizing groupings of length 2.
##                       Check out packages {vioplot} and {see} for alternative ways of plotting split violins/rainclouds

## Warning in x + params$x: longer object length is not a multiple of shorter
## object length
## Warning in x + params$x: longer object length is not a multiple of shorter
## object length
## Warning in x + params$x: longer object length is not a multiple of shorter
## object length

## Warning in compute_layer(..., self = self): Argument 'x' longer than data: some
## values dropped!

ggsave("figures/aim_1_block_synchrony_PFC.png", width = 7, height = 5)

## Warning in setup_data(...): geom_paired_raincloud is only useful for visualizing groupings of length 2.
##                       Check out packages {vioplot} and {see} for alternative ways of plotting split violins/rainclouds

## Warning in x + params$x: longer object length is not a multiple of shorter
## object length
## Warning in x + params$x: longer object length is not a multiple of shorter
## object length
## Warning in x + params$x: longer object length is not a multiple of shorter
## object length

## Warning in compute_layer(..., self = self): Argument 'x' longer than data: some
## values dropped!

DLPFC Graph

synchrony_data %>%
  filter(category == "Real Dyads") %>%
  ggplot(., aes(x = block, y = DLPFC_avg_score, fill = block, col = block)) +
  scale_color_manual(
    name = "Block Type",
    values = c("#C85729", "#127088", "#92874B")
  ) +
  scale_fill_manual(
    name = "Block Type",
    values = c("#C85729", "#127088", "#92874B")
  ) +
  guides(fill = FALSE, col = FALSE) +
  geom_rain(alpha = .5, rain.side = "f2x2", id.long.var = "record_id") +
  ylab("Neural Synchrony Values (DLPFC)") +
  xlab("Play Block") +
  stat_summary(fun.data = "mean_cl_normal", geom = "errorbar", width = 1 / 4, color = "black") +
  stat_summary(fun = "mean", geom = "point", color = "black") +
  theme_apa() +
  theme(text = element_text(size = 15))

## Warning: Duplicated aesthetics after name standardisation: alpha

## Warning in setup_data(...): geom_paired_raincloud is only useful for visualizing groupings of length 2.
##                       Check out packages {vioplot} and {see} for alternative ways of plotting split violins/rainclouds

## Warning in x + params$x: longer object length is not a multiple of shorter
## object length
## Warning in x + params$x: longer object length is not a multiple of shorter
## object length
## Warning in x + params$x: longer object length is not a multiple of shorter
## object length

## Warning in compute_layer(..., self = self): Argument 'x' longer than data: some
## values dropped!

ggsave("figures/aim_1_block_synchrony_DLPFC.png", width = 7, height = 5)

## Warning in setup_data(...): geom_paired_raincloud is only useful for visualizing groupings of length 2.
##                       Check out packages {vioplot} and {see} for alternative ways of plotting split violins/rainclouds

## Warning in x + params$x: longer object length is not a multiple of shorter
## object length
## Warning in x + params$x: longer object length is not a multiple of shorter
## object length
## Warning in x + params$x: longer object length is not a multiple of shorter
## object length

## Warning in compute_layer(..., self = self): Argument 'x' longer than data: some
## values dropped!

VLPFC Graph

synchrony_data %>%
  filter(category == "Real Dyads") %>%
  ggplot(., aes(x = block, y = VLPFC_avg_score, fill = block, col = block)) +
  scale_color_manual(
    name = "Block Type",
    values = c("#C85729", "#127088", "#92874B")
  ) +
  scale_fill_manual(
    name = "Block Type",
    values = c("#C85729", "#127088", "#92874B")
  ) +
  guides(fill = FALSE, col = FALSE) +
  geom_rain(alpha = .5, rain.side = "f2x2", id.long.var = "record_id") +
  ylab("Neural Synchrony Values (VLPFC)") +
  xlab("Play Block") +
  stat_summary(fun.data = "mean_cl_normal", geom = "errorbar", width = 1 / 4, color = "black") +
  stat_summary(fun = "mean", geom = "point", color = "black") +
  theme_apa() +
  theme(text = element_text(size = 15))

## Warning: Duplicated aesthetics after name standardisation: alpha

## Warning in setup_data(...): geom_paired_raincloud is only useful for visualizing groupings of length 2.
##                       Check out packages {vioplot} and {see} for alternative ways of plotting split violins/rainclouds

## Warning in x + params$x: longer object length is not a multiple of shorter
## object length
## Warning in x + params$x: longer object length is not a multiple of shorter
## object length
## Warning in x + params$x: longer object length is not a multiple of shorter
## object length

## Warning in compute_layer(..., self = self): Argument 'x' longer than data: some
## values dropped!

ggsave("figures/aim_1_block_synchrony_VLPFC.png", width = 7, height = 5)

## Warning in setup_data(...): geom_paired_raincloud is only useful for visualizing groupings of length 2.
##                       Check out packages {vioplot} and {see} for alternative ways of plotting split violins/rainclouds

## Warning in x + params$x: longer object length is not a multiple of shorter
## object length
## Warning in x + params$x: longer object length is not a multiple of shorter
## object length
## Warning in x + params$x: longer object length is not a multiple of shorter
## object length

## Warning in compute_layer(..., self = self): Argument 'x' longer than data: some
## values dropped!

Graph per Trial

Block 1: - Trial 1: Markers (fine) - Trial 2: Markers (dots) - Trial 3: Stickers (shapes) - Trial 4: Stickers (strips)

Block 2: - Trial 1: Puzzle set 1 - Trial 2: Puzzle set 2 - Trial 3: Puzzle set 3 - Trial 4: Puzzle set 4

Block 3: - Trial 1: Magnatiles - Trial 2: Human and Animal connectors - Trial 3: Animals, trees, and fences - Trial 4: Cars, trucks

synchrony_data_per_trial$trial <- as.factor(synchrony_data_per_trial$trial)
synchrony_data_per_trial$block <- as.factor(synchrony_data_per_trial$block)

synchrony_data_per_trial$block <- factor(synchrony_data_per_trial$block, levels = c("Block 1", "Block 2", "Block 3"))

synchrony_data_per_trial %>%
  filter(category == "Real Dyads") %>%
  ggplot(., aes(x = trial, y = avg_score, fill = trial, col = trial)) +
  scale_fill_brewer(palette = "Dark2") +
  scale_colour_brewer(palette = "Dark2") +
  guides(fill = FALSE, col = FALSE) +
  geom_rain(alpha = .5, rain.side = "l", id.long.var = "record_id") +
  ylab("Neural Synchrony Values") +
  xlab("Play Block") +
  stat_summary(fun.data = "mean_cl_normal", geom = "errorbar", width = 1 / 4, color = "black", na.rm = TRUE) +
  stat_summary(fun = "mean", geom = "point", color = "black", na.rm = TRUE) +
  theme_apa() +
  theme(text = element_text(size = 20)) +
  facet_wrap(~block) +
  scale_y_continuous(limits = c(0.20, .6), breaks = seq(0.20, .6, by = 0.10)) # Setting Y-axis limits and ticks

## Warning: Duplicated aesthetics after name standardisation: alpha

ggsave("figures/aim_1_block_synchrony_per-trial_PFC.png", width = 14, height = 10)


not_in_stream_df <- synchrony_data_per_trial %>%
  filter(category == "Real Dyads") %>%
  select(record_id, block, trial, DLPFC_avg_score) %>%
  pivot_wider(
    names_from = trial,
    values_from = DLPFC_avg_score
  ) %>%
  rowwise() %>%
  mutate(
    pattern_found = check_sequence(c(`Trial 1`, `Trial 2`, `Trial 3`, `Trial 4`)),
    pattern_found = case_when(
      sum(is.na(c(`Trial 1`, `Trial 2`, `Trial 3`, `Trial 4`))) == 3 ~ TRUE,
      TRUE ~ pattern_found
    )
  ) %>%
  ungroup() %>%
  filter(pattern_found) %>%
  pivot_longer(cols = c("Trial 1":"Trial 4"), values_to = "DLPFC_avg_score", names_to = "trial") %>%
  filter(!is.na(DLPFC_avg_score))


synchrony_data_per_trial %>%
  filter(category == "Real Dyads") %>%
  ggplot(., aes(x = trial, y = DLPFC_avg_score, fill = trial, col = trial)) +
  scale_fill_brewer(palette = "Dark2") +
  scale_colour_brewer(palette = "Dark2") +
  guides(fill = FALSE, col = FALSE) +
  geom_rain(alpha = .5, rain.side = "l", id.long.var = "record_id") +
  geom_point(data = not_in_stream_df, alpha = .5) +
  ylab("Neural Synchrony Values (DLPFC)") +
  xlab("Play Block") +
  stat_summary(fun.data = "mean_cl_normal", geom = "errorbar", width = 1 / 4, color = "black", na.rm = TRUE) +
  stat_summary(fun = "mean", geom = "point", color = "black", na.rm = TRUE) +
  theme_apa() +
  theme(text = element_text(size = 20)) +
  facet_wrap(~block) +
  scale_y_continuous(limits = c(0.20, .6), breaks = seq(0.20, .6, by = 0.10)) # Setting Y-axis limits and ticks

## Warning: Duplicated aesthetics after name standardisation: alpha

## Warning: Removed 51 rows containing non-finite outside the scale range
## (`stat_half_ydensity()`).

## Warning: Removed 51 rows containing non-finite outside the scale range
## (`stat_boxplot()`).

## Warning: Removed 8 rows containing missing values or values outside the scale range
## (`geom_line()`).

## Warning: Removed 51 rows containing missing values or values outside the scale range
## (`geom_point_sorted()`).

ggsave("figures/aim_1_block_synchrony_per-trial_DLPFC.png", width = 14, height = 10)

## Warning: Removed 51 rows containing non-finite outside the scale range
## (`stat_half_ydensity()`).

## Warning: Removed 51 rows containing non-finite outside the scale range
## (`stat_boxplot()`).

## Warning: Removed 8 rows containing missing values or values outside the scale range
## (`geom_line()`).

## Warning: Removed 51 rows containing missing values or values outside the scale range
## (`geom_point_sorted()`).

not_in_stream_df <- synchrony_data_per_trial %>%
  filter(category == "Real Dyads") %>%
  select(record_id, block, trial, VLPFC_avg_score) %>%
  pivot_wider(
    names_from = trial,
    values_from = VLPFC_avg_score
  ) %>%
  rowwise() %>%
  mutate(
    pattern_found = check_sequence(c(`Trial 1`, `Trial 2`, `Trial 3`, `Trial 4`)),
    pattern_found = case_when(
      sum(is.na(c(`Trial 1`, `Trial 2`, `Trial 3`, `Trial 4`))) == 3 ~ TRUE,
      TRUE ~ pattern_found
    )
  ) %>%
  ungroup() %>%
  filter(pattern_found) %>%
  pivot_longer(cols = c("Trial 1":"Trial 4"), values_to = "VLPFC_avg_score", names_to = "trial") %>%
  filter(!is.na(VLPFC_avg_score))

synchrony_data_per_trial %>%
  filter(category == "Real Dyads") %>%
  ggplot(., aes(x = trial, y = VLPFC_avg_score, fill = trial, col = trial)) +
  scale_fill_brewer(palette = "Dark2") +
  scale_colour_brewer(palette = "Dark2") +
  guides(fill = FALSE, col = FALSE) +
  geom_rain(alpha = .5, rain.side = "l", id.long.var = "record_id") +
  geom_point(data = not_in_stream_df, alpha = .5) +
  ylab("Neural Synchrony Values (VLPFC)") +
  xlab("Play Block") +
  stat_summary(fun.data = "mean_cl_normal", geom = "errorbar", width = 1 / 4, color = "black") +
  stat_summary(fun = "mean", geom = "point", color = "black") +
  theme_apa() +
  theme(text = element_text(size = 20)) +
  facet_wrap(~block) +
  scale_y_continuous(limits = c(0.20, .6), breaks = seq(0.20, .6, by = 0.10)) # Setting Y-axis limits and ticks

## Warning: Duplicated aesthetics after name standardisation: alpha

## Warning: Removed 3 rows containing non-finite outside the scale range
## (`stat_half_ydensity()`).

## Warning: Removed 3 rows containing non-finite outside the scale range
## (`stat_boxplot()`).

## Warning: Removed 3 rows containing non-finite outside the scale range
## (`stat_summary()`).
## Removed 3 rows containing non-finite outside the scale range
## (`stat_summary()`).

## Warning: Removed 3 rows containing missing values or values outside the scale range
## (`geom_point_sorted()`).

ggsave("figures/aim_1_block_synchrony_per-trial_VLPFC.png", width = 14, height = 10)

## Warning: Removed 3 rows containing non-finite outside the scale range
## (`stat_half_ydensity()`).

## Warning: Removed 3 rows containing non-finite outside the scale range
## (`stat_boxplot()`).

## Warning: Removed 3 rows containing non-finite outside the scale range
## (`stat_summary()`).
## Removed 3 rows containing non-finite outside the scale range
## (`stat_summary()`).

## Warning: Removed 3 rows containing missing values or values outside the scale range
## (`geom_point_sorted()`).

Aim 2: To examine the relationship between child temperament and dyadic neural synchrony

PFC

Model -block 2 reference

nirs_per_trial_plus_qst$block <- factor(nirs_per_trial_plus_qst$block, levels = c("Block 2", "Block 1", "Block 3"))

model1 <- lmer(
  avg_score ~ c_age_y_centered +
    c_sex + cbq_fear * block + cbq_effortful_control * block + (1 + block | record_id) + (1 | trial),
  data = nirs_per_trial_plus_qst,
  control = lmerControl(
    optimizer = "bobyqa",
    optCtrl = list(maxfun = 5e5)
  )
)

## boundary (singular) fit: see help('isSingular')

summary(model1)

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: 
## avg_score ~ c_age_y_centered + c_sex + cbq_fear * block + cbq_effortful_control *  
##     block + (1 + block | record_id) + (1 | trial)
##    Data: nirs_per_trial_plus_qst
## Control: lmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 5e+05))
## 
## REML criterion at convergence: -4995.8
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.9118 -0.6338 -0.0399  0.6209  3.7242 
## 
## Random effects:
##  Groups    Name         Variance  Std.Dev. Corr       
##  record_id (Intercept)  0.0001712 0.01309             
##            blockBlock 1 0.0002815 0.01678  -1.00      
##            blockBlock 3 0.0001176 0.01084  -1.00  1.00
##  trial     (Intercept)  0.0000000 0.00000             
##  Residual               0.0011164 0.03341             
## Number of obs: 1303, groups:  record_id, 109; trial, 4
## 
## Fixed effects:
##                                      Estimate Std. Error         df t value
## (Intercept)                         3.348e-01  1.845e-02  1.115e+02  18.141
## c_age_y_centered                    4.687e-04  8.481e-04  3.245e+02   0.553
## c_sexF                              4.212e-04  2.070e-03  3.257e+02   0.203
## cbq_fear                           -2.019e-04  1.934e-03  1.101e+02  -0.104
## blockBlock 1                        1.937e-02  2.508e-02  1.154e+02   0.772
## blockBlock 3                        2.933e-02  2.249e-02  1.597e+02   1.304
## cbq_effortful_control               5.746e-03  3.258e-03  1.135e+02   1.763
## cbq_fear:blockBlock 1              -2.632e-04  2.639e-03  1.152e+02  -0.100
## cbq_fear:blockBlock 3              -1.205e-04  2.376e-03  1.621e+02  -0.051
## blockBlock 1:cbq_effortful_control -7.652e-03  4.406e-03  1.156e+02  -1.737
## blockBlock 3:cbq_effortful_control -9.727e-03  3.951e-03  1.600e+02  -2.462
##                                    Pr(>|t|)    
## (Intercept)                          <2e-16 ***
## c_age_y_centered                     0.5809    
## c_sexF                               0.8389    
## cbq_fear                             0.9171    
## blockBlock 1                         0.4415    
## blockBlock 3                         0.1941    
## cbq_effortful_control                0.0805 .  
## cbq_fear:blockBlock 1                0.9207    
## cbq_fear:blockBlock 3                0.9596    
## blockBlock 1:cbq_effortful_control   0.0851 .  
## blockBlock 3:cbq_effortful_control   0.0149 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) c_g_y_ c_sexF cbq_fr blckB1 blckB3 cbq_f_ cb_:B1 cb_:B3
## c_g_y_cntrd -0.084                                                        
## c_sexF       0.073 -0.252                                                 
## cbq_fear    -0.403  0.042 -0.020                                          
## blockBlock1 -0.807  0.000 -0.001  0.325                                   
## blockBlock3 -0.759  0.000 -0.001  0.306  0.612                            
## cbq_ffrtfl_ -0.902  0.088 -0.132 -0.012  0.724  0.681                     
## cbq_fr:blB1  0.324  0.000  0.000 -0.810 -0.401 -0.245  0.014              
## cbq_fr:blB3  0.304  0.001  0.004 -0.760 -0.245 -0.401  0.013  0.610       
## blckB1:cb__  0.728  0.000  0.001  0.014 -0.902 -0.552 -0.802 -0.017 -0.011
## blckB3:cb__  0.685  0.000 -0.001  0.013 -0.552 -0.901 -0.755 -0.011 -0.020
##             bB1:__
## c_g_y_cntrd       
## c_sexF            
## cbq_fear          
## blockBlock1       
## blockBlock3       
## cbq_ffrtfl_       
## cbq_fr:blB1       
## cbq_fr:blB3       
## blckB1:cb__       
## blckB3:cb__  0.611
## optimizer (bobyqa) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')

# emm = emmeans(model1, ~ cbq_effortful_control*block )
# pairs(emm)
# # or for simple comparisons
# #pairs(emm, simple = "each")

anova(model1, type = "III")

report(model1)

## boundary (singular) fit: see help('isSingular')

## Random effect variances not available. Returned R2 does not account for random effects.

## boundary (singular) fit: see help('isSingular')

## Random effect variances not available. Returned R2 does not account for random effects.

## We fitted a linear mixed model (estimated using REML and BOBYQA optimizer) to
## predict avg_score with c_age_y_centered, c_sex, cbq_fear, block and
## cbq_effortful_control (formula: avg_score ~ c_age_y_centered + c_sex + cbq_fear
## * block + cbq_effortful_control * block). The model included block as random
## effects (formula: list(~1 + block | record_id, ~1 | trial)). The model's
## explanatory power related to the fixed effects alone (marginal R2) is 0.09. The
## model's intercept, corresponding to c_age_y_centered = 0, c_sex = M, cbq_fear =
## 0, block = Block 2 and cbq_effortful_control = 0, is at 0.33 (95% CI [0.30,
## 0.37], t(1284) = 18.14, p < .001). Within this model:
## 
##   - The effect of c age y centered is statistically non-significant and positive
## (beta = 4.69e-04, 95% CI [-1.20e-03, 2.13e-03], t(1284) = 0.55, p = 0.581; Std.
## beta = 0.02, 95% CI [-0.04, 0.07])
##   - The effect of c sex [F] is statistically non-significant and positive (beta =
## 4.21e-04, 95% CI [-3.64e-03, 4.48e-03], t(1284) = 0.20, p = 0.839; Std. beta =
## 0.01, 95% CI [-0.10, 0.13])
##   - The effect of cbq fear is statistically non-significant and negative (beta =
## -2.02e-04, 95% CI [-4.00e-03, 3.59e-03], t(1284) = -0.10, p = 0.917; Std. beta
## = -5.94e-03, 95% CI [-0.12, 0.11])
##   - The effect of block [Block 1] is statistically non-significant and positive
## (beta = 0.02, 95% CI [-0.03, 0.07], t(1284) = 0.77, p = 0.440; Std. beta =
## -0.60, 95% CI [-0.75, -0.44])
##   - The effect of block [Block 3] is statistically non-significant and positive
## (beta = 0.03, 95% CI [-0.01, 0.07], t(1284) = 1.30, p = 0.192; Std. beta =
## -0.60, 95% CI [-0.74, -0.46])
##   - The effect of cbq effortful control is statistically non-significant and
## positive (beta = 5.75e-03, 95% CI [-6.46e-04, 0.01], t(1284) = 1.76, p = 0.078;
## Std. beta = 0.10, 95% CI [-0.01, 0.21])
##   - The effect of cbq fear × block [Block 1] is statistically non-significant and
## negative (beta = -2.63e-04, 95% CI [-5.44e-03, 4.91e-03], t(1284) = -0.10, p =
## 0.921; Std. beta = -7.74e-03, 95% CI [-0.16, 0.14])
##   - The effect of cbq fear × block [Block 3] is statistically non-significant and
## negative (beta = -1.20e-04, 95% CI [-4.78e-03, 4.54e-03], t(1284) = -0.05, p =
## 0.960; Std. beta = -3.54e-03, 95% CI [-0.14, 0.13])
##   - The effect of block [Block 1] × cbq effortful control is statistically
## non-significant and negative (beta = -7.65e-03, 95% CI [-0.02, 9.92e-04],
## t(1284) = -1.74, p = 0.083; Std. beta = -0.14, 95% CI [-0.29, 0.02])
##   - The effect of block [Block 3] × cbq effortful control is statistically
## significant and negative (beta = -9.73e-03, 95% CI [-0.02, -1.98e-03], t(1284)
## = -2.46, p = 0.014; Std. beta = -0.17, 95% CI [-0.31, -0.03])
## 
## Standardized parameters were obtained by fitting the model on a standardized
## version of the dataset. 95% Confidence Intervals (CIs) and p-values were
## computed using a Wald t-distribution approximation.

performance::check_outliers(model1)

## OK: No outliers detected.
## - Based on the following method and threshold: cook (0.9).
## - For variable: (Whole model)

nirs_per_trial_plus_qst$block_numeric <- as.numeric(nirs_per_trial_plus_qst$block)

model1b <- lmer(
  avg_score ~ c_age_y_centered +
    c_sex + cbq_fear * block_numeric + cbq_effortful_control * block_numeric + (1 + block | record_id) + (1 | trial),
  data = nirs_per_trial_plus_qst,
  control = lmerControl(
    optimizer = "bobyqa",
    optCtrl = list(maxfun = 5e5)
  )
)

## boundary (singular) fit: see help('isSingular')

summary(model1b)

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: avg_score ~ c_age_y_centered + c_sex + cbq_fear * block_numeric +  
##     cbq_effortful_control * block_numeric + (1 + block | record_id) +  
##     (1 | trial)
##    Data: nirs_per_trial_plus_qst
## Control: lmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 5e+05))
## 
## REML criterion at convergence: -5000.6
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.9192 -0.6381 -0.0496  0.6412  3.7989 
## 
## Random effects:
##  Groups    Name         Variance  Std.Dev. Corr       
##  record_id (Intercept)  0.0002026 0.01424             
##            blockBlock 1 0.0003939 0.01985  -1.00      
##            blockBlock 3 0.0001195 0.01093  -1.00  1.00
##  trial     (Intercept)  0.0000000 0.00000             
##  Residual               0.0011220 0.03350             
## Number of obs: 1303, groups:  record_id, 109; trial, 4
## 
## Fixed effects:
##                                       Estimate Std. Error         df t value
## (Intercept)                          3.232e-01  2.460e-02  2.249e+02  13.138
## c_age_y_centered                     4.416e-04  8.475e-04  3.692e+02   0.521
## c_sexF                               4.078e-04  2.069e-03  3.710e+02   0.197
## cbq_fear                            -3.142e-04  2.588e-03  2.239e+02  -0.121
## block_numeric                        1.407e-02  1.096e-02  2.720e+02   1.284
## cbq_effortful_control                8.780e-03  4.334e-03  2.269e+02   2.026
## cbq_fear:block_numeric              -1.231e-05  1.158e-03  2.741e+02  -0.011
## block_numeric:cbq_effortful_control -4.534e-03  1.925e-03  2.721e+02  -2.356
##                                     Pr(>|t|)    
## (Intercept)                           <2e-16 ***
## c_age_y_centered                      0.6026    
## c_sexF                                0.8438    
## cbq_fear                              0.9035    
## block_numeric                         0.2002    
## cbq_effortful_control                 0.0439 *  
## cbq_fear:block_numeric                0.9915    
## block_numeric:cbq_effortful_control   0.0192 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) c_g_y_ c_sexF cbq_fr blck_n cbq_f_ cbq_:_
## c_g_y_cntrd -0.063                                          
## c_sexF       0.054 -0.252                                   
## cbq_fear    -0.402  0.031 -0.018                            
## block_numrc -0.931  0.000 -0.001  0.375                     
## cbq_ffrtfl_ -0.902  0.066 -0.097 -0.016  0.837              
## cbq_fr:blc_  0.373  0.001  0.005 -0.933 -0.401  0.018       
## blck_nmr:__  0.839  0.000 -0.002  0.018 -0.901 -0.928 -0.021
## optimizer (bobyqa) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')

report(model1b)

## boundary (singular) fit: see help('isSingular')

## Random effect variances not available. Returned R2 does not account for random effects.

## boundary (singular) fit: see help('isSingular')

## Random effect variances not available. Returned R2 does not account for random effects.

## We fitted a linear mixed model (estimated using REML and BOBYQA optimizer) to
## predict avg_score with c_age_y_centered, c_sex, cbq_fear, block_numeric and
## cbq_effortful_control (formula: avg_score ~ c_age_y_centered + c_sex + cbq_fear
## * block_numeric + cbq_effortful_control * block_numeric). The model included
## block as random effects (formula: list(~1 + block | record_id, ~1 | trial)).
## The model's explanatory power related to the fixed effects alone (marginal R2)
## is 0.06. The model's intercept, corresponding to c_age_y_centered = 0, c_sex =
## M, cbq_fear = 0, block_numeric = 0 and cbq_effortful_control = 0, is at 0.32
## (95% CI [0.27, 0.37], t(1287) = 13.14, p < .001). Within this model:
## 
##   - The effect of c age y centered is statistically non-significant and positive
## (beta = 4.42e-04, 95% CI [-1.22e-03, 2.10e-03], t(1287) = 0.52, p = 0.602; Std.
## beta = 0.01, 95% CI [-0.04, 0.07])
##   - The effect of c sex [F] is statistically non-significant and positive (beta =
## 4.08e-04, 95% CI [-3.65e-03, 4.47e-03], t(1287) = 0.20, p = 0.844; Std. beta =
## 0.01, 95% CI [-0.10, 0.12])
##   - The effect of cbq fear is statistically non-significant and negative (beta =
## -3.14e-04, 95% CI [-5.39e-03, 4.76e-03], t(1287) = -0.12, p = 0.903; Std. beta
## = -9.96e-03, 95% CI [-0.06, 0.04])
##   - The effect of block numeric is statistically non-significant and positive
## (beta = 0.01, 95% CI [-7.42e-03, 0.04], t(1287) = 1.28, p = 0.199; Std. beta =
## -0.22, 95% CI [-0.27, -0.16])
##   - The effect of cbq effortful control is statistically significant and positive
## (beta = 8.78e-03, 95% CI [2.78e-04, 0.02], t(1287) = 2.03, p = 0.043; Std. beta
## = -4.90e-03, 95% CI [-0.06, 0.05])
##   - The effect of cbq fear × block numeric is statistically non-significant and
## negative (beta = -1.23e-05, 95% CI [-2.28e-03, 2.26e-03], t(1287) = -0.01, p =
## 0.992; Std. beta = -2.96e-04, 95% CI [-0.05, 0.05])
##   - The effect of block numeric × cbq effortful control is statistically
## significant and negative (beta = -4.53e-03, 95% CI [-8.31e-03, -7.58e-04],
## t(1287) = -2.36, p = 0.019; Std. beta = -0.07, 95% CI [-0.12, -0.01])
## 
## Standardized parameters were obtained by fitting the model on a standardized
## version of the dataset. 95% Confidence Intervals (CIs) and p-values were
## computed using a Wald t-distribution approximation.

Report

tab_model(model1, p.val = "kr", show.df = TRUE)

	avg score
Predictors	Estimates	CI	p	df
(Intercept)	0.33	0.30 – 0.37	<0.001	106.76
c age y centered	0.00	-0.00 – 0.00	0.589	103.81
c sex [F]	0.00	-0.00 – 0.00	0.842	103.82
cbq fear	-0.00	-0.00 – 0.00	0.917	105.48
block [Block 1]	0.02	-0.03 – 0.07	0.442	105.86
block [Block 3]	0.03	-0.02 – 0.07	0.195	105.51
cbq effortful control	0.01	-0.00 – 0.01	0.081	108.57
cbq fear × block [Block 1]	-0.00	-0.01 – 0.00	0.921	105.69
cbq fear × block [Block 3]	-0.00	-0.00 – 0.00	0.960	106.51
block [Block 1] × cbq effortful control	-0.01	-0.02 – 0.00	0.085	105.96
block [Block 3] × cbq effortful control	-0.01	-0.02 – -0.00	0.015	105.60
Random Effects
σ2	0.00
τ00 record_id	0.00
τ00 trial	0.00
τ11 record_id.blockBlock 1	0.00
τ11 record_id.blockBlock 3	0.00
ρ01 record_id.blockBlock 1	-1.00
ρ01 record_id.blockBlock 3	-1.00
N record_id	109
N trial	4
Observations	1303
Marginal R2 / Conditional R2	0.089 / NA

Simple intercepts and slopes

sim_slopes(model1,
  pred = cbq_effortful_control, modx = block, cond.int = TRUE,
  centered = "none", johnson_neyman = FALSE
)

## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')

## Warning: Model failed to converge with 1 negative eigenvalue: -4.1e+00

## SIMPLE SLOPES ANALYSIS
## 
## When block = Block 2: 
## 
##                                        Est.   S.E.   t val.      p
## ------------------------------------ ------ ------ -------- ------
## Slope of cbq_effortful_control         0.01   0.00     1.76   0.08
## Conditional intercept                  0.33   0.02    18.14   0.00
## 
## When block = Block 1: 
## 
##                                         Est.   S.E.   t val.      p
## ------------------------------------ ------- ------ -------- ------
## Slope of cbq_effortful_control         -0.00   0.00    -0.72   0.47
## Conditional intercept                   0.35   0.01    23.71   0.00
## 
## When block = Block 3: 
## 
##                                         Est.   S.E.   t val.      p
## ------------------------------------ ------- ------ -------- ------
## Slope of cbq_effortful_control         -0.00   0.00    -1.53   0.13
## Conditional intercept                   0.36   0.01    24.76   0.00

Interaction Plot

interact_plot(model1,
  pred = cbq_effortful_control, modx = block,
  x.label = "Effortful Control", y.label = "PFC Synchrony Scores",
  legend.main = "Block"
)

Plot for interaction: EC*Block

interaction_eff1 <- ggpredict(model1, c("cbq_effortful_control", "block"))
interaction_eff1

ggplot(interaction_eff1, aes(x, predicted, color = group)) +
  geom_line(linewidth = 1.2) +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high, fill = group), alpha = 0.2, color = NA) +
  theme_classic() +
  scale_color_manual(
    name = "Block Type",
    values = natparks.pals("Saguaro", 3),
    labels = c(
      "Block 1" = "Baseline",
      "Block 2" = "Stressful",
      "Block 3" = "Recovery"
    )
  ) +
  scale_fill_manual(
    name = "Block Type",
    values = natparks.pals("Saguaro", 3),
    labels = c(
      "Block 1" = "Baseline",
      "Block 2" = "Stressful",
      "Block 3" = "Recovery"
    )
  ) +
  scale_x_continuous(
    name = "Effortful Control (CBQ scores)"
  ) +
  theme(
    axis.text = element_text(size = 19),
    axis.title = element_text(size = 20),
    legend.text = element_text(size = 19),
    legend.title = element_text(face = "bold", size = 19),
    legend.title.align = 0.5
  ) +
  scale_y_continuous(name = "PFC Neural Synchrony (predicted)", limits = c(0.29, .42), breaks = seq(0.29, .42, by = 0.02)) # Setting Y-axis limits and ticks

## Warning: The `legend.title.align` argument of `theme()` is deprecated as of ggplot2
## 3.5.0.
## ℹ Please use theme(legend.title = element_text(hjust)) instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

ggsave("figures/aim_2_EC-Block_PFC.png", width = 10, height = 8)

Plot for interaction: F*Block

interaction_eff1b <- ggpredict(model1, c("cbq_fear", "block"))
interaction_eff1b

ggplot(interaction_eff1b, aes(x, predicted, color = group)) +
  geom_line(linewidth = 1.2) +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high, fill = group), alpha = 0.2, color = NA) +
  theme_classic() +
  scale_color_manual(
    name = "Block Type",
    values = natparks.pals("Saguaro", 3),
    labels = c(
      "Block 1" = "Baseline",
      "Block 2" = "Stressful",
      "Block 3" = "Recovery"
    )
  ) +
  scale_fill_manual(
    name = "Block Type",
    values = natparks.pals("Saguaro", 3),
    labels = c(
      "Block 1" = "Baseline",
      "Block 2" = "Stressful",
      "Block 3" = "Recovery"
    )
  ) +
  scale_x_continuous(
    name = "Fear (CBQ scores)"
  ) +
  theme(
    axis.text = element_text(size = 19),
    axis.title = element_text(size = 20),
    legend.text = element_text(size = 19),
    legend.title = element_text(face = "bold", size = 19),
    legend.title.align = 0.5
  ) +
  scale_y_continuous(name = "PFC Neural Synchrony (predicted)", limits = c(0.29, .42), breaks = seq(0.29, .42, by = 0.02)) # Setting Y-axis limits and ticks

ggsave("figures/aim_2_F-Block_PFC.png", width = 10, height = 8)

DLPFC

Model-block 2 reference

nirs_per_trial_plus_qst$block <- factor(nirs_per_trial_plus_qst$block, levels = c("Block 2", "Block 1", "Block 3"))

model2 <- lmer(
  DLPFC_avg_score ~ c_age_y_centered +
    c_sex + cbq_fear * block + cbq_effortful_control * block + (1 + block | record_id) + (1 | trial),
  data = nirs_per_trial_plus_qst,
  control = lmerControl(
    optimizer = "bobyqa",
    optCtrl = list(maxfun = 5e5)
  )
)

## boundary (singular) fit: see help('isSingular')

summary(model2)

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: DLPFC_avg_score ~ c_age_y_centered + c_sex + cbq_fear * block +  
##     cbq_effortful_control * block + (1 + block | record_id) +      (1 | trial)
##    Data: nirs_per_trial_plus_qst
## Control: lmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 5e+05))
## 
## REML criterion at convergence: -3666.8
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.2707 -0.6877 -0.0141  0.6253  3.6541 
## 
## Random effects:
##  Groups    Name         Variance  Std.Dev. Corr       
##  record_id (Intercept)  0.0003342 0.01828             
##            blockBlock 1 0.0003763 0.01940  -1.00      
##            blockBlock 3 0.0005124 0.02264  -1.00  1.00
##  trial     (Intercept)  0.0000000 0.00000             
##  Residual               0.0028132 0.05304             
## Number of obs: 1254, groups:  record_id, 109; trial, 4
## 
## Fixed effects:
##                                      Estimate Std. Error         df t value
## (Intercept)                         3.139e-01  2.831e-02  1.030e+02  11.087
## c_age_y_centered                    2.693e-03  1.347e-03  4.518e+02   2.000
## c_sexF                             -5.927e-04  3.290e-03  4.525e+02  -0.180
## cbq_fear                            7.029e-03  2.985e-03  1.027e+02   2.355
## blockBlock 1                        7.367e-03  3.701e-02  1.296e+02   0.199
## blockBlock 3                        8.359e-02  3.827e-02  1.123e+02   2.184
## cbq_effortful_control               4.730e-03  4.996e-03  1.048e+02   0.947
## cbq_fear:blockBlock 1              -5.165e-03  3.906e-03  1.302e+02  -1.322
## cbq_fear:blockBlock 3              -9.675e-03  4.060e-03  1.151e+02  -2.383
## blockBlock 1:cbq_effortful_control -1.522e-03  6.489e-03  1.292e+02  -0.235
## blockBlock 3:cbq_effortful_control -1.346e-02  6.734e-03  1.132e+02  -1.999
##                                    Pr(>|t|)    
## (Intercept)                          <2e-16 ***
## c_age_y_centered                     0.0461 *  
## c_sexF                               0.8571    
## cbq_fear                             0.0204 *  
## blockBlock 1                         0.8425    
## blockBlock 3                         0.0310 *  
## cbq_effortful_control                0.3459    
## cbq_fear:blockBlock 1                0.1884    
## cbq_fear:blockBlock 3                0.0188 *  
## blockBlock 1:cbq_effortful_control   0.8149    
## blockBlock 3:cbq_effortful_control   0.0480 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) c_g_y_ c_sexF cbq_fr blckB1 blckB3 cbq_f_ cb_:B1 cb_:B3
## c_g_y_cntrd -0.081                                                        
## c_sexF       0.071 -0.281                                                 
## cbq_fear    -0.407  0.036 -0.018                                          
## blockBlock1 -0.772 -0.008  0.003  0.315                                   
## blockBlock3 -0.788 -0.006  0.001  0.321  0.616                            
## cbq_ffrtfl_ -0.900  0.090 -0.133 -0.011  0.691  0.705                     
## cbq_fr:blB1  0.314  0.009 -0.003 -0.777 -0.407 -0.251  0.013              
## cbq_fr:blB3  0.318  0.008  0.001 -0.788 -0.249 -0.399  0.012  0.616       
## blckB1:cb__  0.696  0.005 -0.001  0.012 -0.902 -0.556 -0.768 -0.012 -0.009
## blckB3:cb__  0.708  0.003 -0.001  0.012 -0.554 -0.900 -0.781 -0.009 -0.024
##             bB1:__
## c_g_y_cntrd       
## c_sexF            
## cbq_fear          
## blockBlock1       
## blockBlock3       
## cbq_ffrtfl_       
## cbq_fr:blB1       
## cbq_fr:blB3       
## blckB1:cb__       
## blckB3:cb__  0.615
## optimizer (bobyqa) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')

report(model2)

## boundary (singular) fit: see help('isSingular')

## Random effect variances not available. Returned R2 does not account for random effects.

## boundary (singular) fit: see help('isSingular')

## Random effect variances not available. Returned R2 does not account for random effects.

## We fitted a linear mixed model (estimated using REML and BOBYQA optimizer) to
## predict DLPFC_avg_score with c_age_y_centered, c_sex, cbq_fear, block and
## cbq_effortful_control (formula: DLPFC_avg_score ~ c_age_y_centered + c_sex +
## cbq_fear * block + cbq_effortful_control * block). The model included block as
## random effects (formula: list(~1 + block | record_id, ~1 | trial)). The model's
## explanatory power related to the fixed effects alone (marginal R2) is 0.05. The
## model's intercept, corresponding to c_age_y_centered = 0, c_sex = M, cbq_fear =
## 0, block = Block 2 and cbq_effortful_control = 0, is at 0.31 (95% CI [0.26,
## 0.37], t(1235) = 11.09, p < .001). Within this model:
## 
##   - The effect of c age y centered is statistically significant and positive
## (beta = 2.69e-03, 95% CI [5.08e-05, 5.34e-03], t(1235) = 2.00, p = 0.046; Std.
## beta = 0.06, 95% CI [1.10e-03, 0.12])
##   - The effect of c sex [F] is statistically non-significant and negative (beta =
## -5.93e-04, 95% CI [-7.05e-03, 5.86e-03], t(1235) = -0.18, p = 0.857; Std. beta
## = -0.01, 95% CI [-0.13, 0.11])
##   - The effect of cbq fear is statistically significant and positive (beta =
## 7.03e-03, 95% CI [1.17e-03, 0.01], t(1235) = 2.35, p = 0.019; Std. beta = 0.13,
## 95% CI [0.02, 0.25])
##   - The effect of block [Block 1] is statistically non-significant and positive
## (beta = 7.37e-03, 95% CI [-0.07, 0.08], t(1235) = 0.20, p = 0.842; Std. beta =
## -0.38, 95% CI [-0.52, -0.23])
##   - The effect of block [Block 3] is statistically significant and positive (beta
## = 0.08, 95% CI [8.50e-03, 0.16], t(1235) = 2.18, p = 0.029; Std. beta = -0.44,
## 95% CI [-0.59, -0.29])
##   - The effect of cbq effortful control is statistically non-significant and
## positive (beta = 4.73e-03, 95% CI [-5.07e-03, 0.01], t(1235) = 0.95, p = 0.344;
## Std. beta = 0.05, 95% CI [-0.06, 0.17])
##   - The effect of cbq fear × block [Block 1] is statistically non-significant and
## negative (beta = -5.17e-03, 95% CI [-0.01, 2.50e-03], t(1235) = -1.32, p =
## 0.186; Std. beta = -0.10, 95% CI [-0.24, 0.05])
##   - The effect of cbq fear × block [Block 3] is statistically significant and
## negative (beta = -9.67e-03, 95% CI [-0.02, -1.71e-03], t(1235) = -2.38, p =
## 0.017; Std. beta = -0.18, 95% CI [-0.34, -0.03])
##   - The effect of block [Block 1] × cbq effortful control is statistically
## non-significant and negative (beta = -1.52e-03, 95% CI [-0.01, 0.01], t(1235) =
## -0.23, p = 0.815; Std. beta = -0.02, 95% CI [-0.16, 0.13])
##   - The effect of block [Block 3] × cbq effortful control is statistically
## significant and negative (beta = -0.01, 95% CI [-0.03, -2.52e-04], t(1235) =
## -2.00, p = 0.046; Std. beta = -0.15, 95% CI [-0.31, -2.89e-03])
## 
## Standardized parameters were obtained by fitting the model on a standardized
## version of the dataset. 95% Confidence Intervals (CIs) and p-values were
## computed using a Wald t-distribution approximation.

emm <- emmeans(model2, ~ cbq_effortful_control * block)
pairs(emm)

##  contrast                                                                                     
##  cbq_effortful_control5.17144338118022 Block 2 - cbq_effortful_control5.17144338118022 Block 1
##  cbq_effortful_control5.17144338118022 Block 2 - cbq_effortful_control5.17144338118022 Block 3
##  cbq_effortful_control5.17144338118022 Block 1 - cbq_effortful_control5.17144338118022 Block 3
##  estimate      SE  df t.ratio p.value
##   0.02097 0.00413 105   5.079  <.0001
##   0.02438 0.00427 105   5.714  <.0001
##   0.00341 0.00368 105   0.926  0.6254
## 
## Results are averaged over the levels of: c_sex 
## Degrees-of-freedom method: kenward-roger 
## P value adjustment: tukey method for comparing a family of 3 estimates

# or for simple comparisons
# pairs(emm, simple = "each")

Model-block 3 reference

nirs_per_trial_plus_qst$block <- factor(nirs_per_trial_plus_qst$block, levels = c("Block 3", "Block 1", "Block 2"))

model2b <- lmer(
  DLPFC_avg_score ~ c_age_y_centered +
    c_sex + cbq_fear * block + cbq_effortful_control * block + (1 + block | record_id) + (1 | trial),
  data = nirs_per_trial_plus_qst,
  control = lmerControl(
    optimizer = "bobyqa",
    optCtrl = list(maxfun = 5e5)
  )
)

## boundary (singular) fit: see help('isSingular')

summary(model2b)

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: DLPFC_avg_score ~ c_age_y_centered + c_sex + cbq_fear * block +  
##     cbq_effortful_control * block + (1 + block | record_id) +      (1 | trial)
##    Data: nirs_per_trial_plus_qst
## Control: lmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 5e+05))
## 
## REML criterion at convergence: -3666.8
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.2707 -0.6877 -0.0141  0.6253  3.6541 
## 
## Random effects:
##  Groups    Name         Variance  Std.Dev. Corr       
##  record_id (Intercept)  1.898e-05 0.004357            
##            blockBlock 1 1.049e-05 0.003239 -1.00      
##            blockBlock 2 5.124e-04 0.022637 -1.00  1.00
##  trial     (Intercept)  0.000e+00 0.000000            
##  Residual               2.813e-03 0.053039            
## Number of obs: 1254, groups:  record_id, 109; trial, 4
## 
## Fixed effects:
##                                      Estimate Std. Error         df t value
## (Intercept)                         3.975e-01  2.365e-02  5.378e+02  16.806
## c_age_y_centered                    2.693e-03  1.347e-03  4.518e+02   2.000
## c_sexF                             -5.927e-04  3.290e-03  4.525e+02  -0.180
## cbq_fear                           -2.646e-03  2.507e-03  5.337e+02  -1.056
## blockBlock 1                       -7.622e-02  3.299e-02  8.443e+02  -2.310
## blockBlock 2                       -8.359e-02  3.827e-02  1.123e+02  -2.184
## cbq_effortful_control              -8.734e-03  4.215e-03  5.632e+02  -2.072
## cbq_fear:blockBlock 1               4.510e-03  3.493e-03  8.469e+02   1.291
## cbq_fear:blockBlock 2               9.675e-03  4.060e-03  1.151e+02   2.383
## blockBlock 1:cbq_effortful_control  1.194e-02  5.807e-03  8.457e+02   2.056
## blockBlock 2:cbq_effortful_control  1.346e-02  6.734e-03  1.132e+02   1.999
##                                    Pr(>|t|)    
## (Intercept)                          <2e-16 ***
## c_age_y_centered                     0.0461 *  
## c_sexF                               0.8571    
## cbq_fear                             0.2916    
## blockBlock 1                         0.0211 *  
## blockBlock 2                         0.0310 *  
## cbq_effortful_control                0.0387 *  
## cbq_fear:blockBlock 1                0.1971    
## cbq_fear:blockBlock 2                0.0188 *  
## blockBlock 1:cbq_effortful_control   0.0400 *  
## blockBlock 2:cbq_effortful_control   0.0480 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) c_g_y_ c_sexF cbq_fr blckB1 blckB2 cbq_f_ cb_:B1 cb_:B2
## c_g_y_cntrd -0.106                                                        
## c_sexF       0.085 -0.281                                                 
## cbq_fear    -0.392  0.055 -0.020                                          
## blockBlock1 -0.701 -0.002  0.002  0.274                                   
## blockBlock2 -0.675  0.006 -0.001  0.265  0.469                            
## cbq_ffrtfl_ -0.898  0.111 -0.160 -0.030  0.625  0.602                     
## cbq_fr:blB1  0.274  0.001 -0.004 -0.711 -0.399 -0.184  0.027              
## cbq_fr:blB2  0.265 -0.008 -0.001 -0.681 -0.183 -0.399  0.024  0.473       
## blckB1:cb__  0.633  0.002  0.000  0.027 -0.901 -0.423 -0.699 -0.023 -0.017
## blckB2:cb__  0.610 -0.003  0.001  0.024 -0.423 -0.900 -0.672 -0.017 -0.024
##             bB1:__
## c_g_y_cntrd       
## c_sexF            
## cbq_fear          
## blockBlock1       
## blockBlock2       
## cbq_ffrtfl_       
## cbq_fr:blB1       
## cbq_fr:blB2       
## blckB1:cb__       
## blckB2:cb__  0.473
## optimizer (bobyqa) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')

report(model2b)

## boundary (singular) fit: see help('isSingular')

## Random effect variances not available. Returned R2 does not account for random effects.

## boundary (singular) fit: see help('isSingular')

## Random effect variances not available. Returned R2 does not account for random effects.

## We fitted a linear mixed model (estimated using REML and BOBYQA optimizer) to
## predict DLPFC_avg_score with c_age_y_centered, c_sex, cbq_fear, block and
## cbq_effortful_control (formula: DLPFC_avg_score ~ c_age_y_centered + c_sex +
## cbq_fear * block + cbq_effortful_control * block). The model included block as
## random effects (formula: list(~1 + block | record_id, ~1 | trial)). The model's
## explanatory power related to the fixed effects alone (marginal R2) is 0.05. The
## model's intercept, corresponding to c_age_y_centered = 0, c_sex = M, cbq_fear =
## 0, block = Block 3 and cbq_effortful_control = 0, is at 0.40 (95% CI [0.35,
## 0.44], t(1235) = 16.81, p < .001). Within this model:
## 
##   - The effect of c age y centered is statistically significant and positive
## (beta = 2.69e-03, 95% CI [5.08e-05, 5.34e-03], t(1235) = 2.00, p = 0.046; Std.
## beta = 0.06, 95% CI [1.10e-03, 0.12])
##   - The effect of c sex [F] is statistically non-significant and negative (beta =
## -5.93e-04, 95% CI [-7.05e-03, 5.86e-03], t(1235) = -0.18, p = 0.857; Std. beta
## = -0.01, 95% CI [-0.13, 0.11])
##   - The effect of cbq fear is statistically non-significant and negative (beta =
## -2.65e-03, 95% CI [-7.56e-03, 2.27e-03], t(1235) = -1.06, p = 0.291; Std. beta
## = -0.05, 95% CI [-0.14, 0.04])
##   - The effect of block [Block 1] is statistically significant and negative (beta
## = -0.08, 95% CI [-0.14, -0.01], t(1235) = -2.31, p = 0.021; Std. beta = 0.06,
## 95% CI [-0.07, 0.19])
##   - The effect of block [Block 2] is statistically significant and negative (beta
## = -0.08, 95% CI [-0.16, -8.50e-03], t(1235) = -2.18, p = 0.029; Std. beta =
## 0.44, 95% CI [0.29, 0.59])
##   - The effect of cbq effortful control is statistically significant and negative
## (beta = -8.73e-03, 95% CI [-0.02, -4.65e-04], t(1235) = -2.07, p = 0.038; Std.
## beta = -0.10, 95% CI [-0.19, -5.33e-03])
##   - The effect of cbq fear × block [Block 1] is statistically non-significant and
## positive (beta = 4.51e-03, 95% CI [-2.34e-03, 0.01], t(1235) = 1.29, p = 0.197;
## Std. beta = 0.09, 95% CI [-0.04, 0.22])
##   - The effect of cbq fear × block [Block 2] is statistically significant and
## positive (beta = 9.67e-03, 95% CI [1.71e-03, 0.02], t(1235) = 2.38, p = 0.017;
## Std. beta = 0.18, 95% CI [0.03, 0.34])
##   - The effect of block [Block 1] × cbq effortful control is statistically
## significant and positive (beta = 0.01, 95% CI [5.49e-04, 0.02], t(1235) = 2.06,
## p = 0.040; Std. beta = 0.14, 95% CI [6.29e-03, 0.27])
##   - The effect of block [Block 2] × cbq effortful control is statistically
## significant and positive (beta = 0.01, 95% CI [2.52e-04, 0.03], t(1235) = 2.00,
## p = 0.046; Std. beta = 0.15, 95% CI [2.89e-03, 0.31])
## 
## Standardized parameters were obtained by fitting the model on a standardized
## version of the dataset. 95% Confidence Intervals (CIs) and p-values were
## computed using a Wald t-distribution approximation.

emm <- emmeans(model2b, ~ cbq_effortful_control * block)
pairs(emm)

##  contrast                                                                                     
##  cbq_effortful_control5.17144338118022 Block 3 - cbq_effortful_control5.17144338118022 Block 1
##  cbq_effortful_control5.17144338118022 Block 3 - cbq_effortful_control5.17144338118022 Block 2
##  cbq_effortful_control5.17144338118022 Block 1 - cbq_effortful_control5.17144338118022 Block 2
##  estimate      SE  df t.ratio p.value
##  -0.00341 0.00368 105  -0.926  0.6254
##  -0.02438 0.00427 105  -5.714  <.0001
##  -0.02097 0.00413 105  -5.079  <.0001
## 
## Results are averaged over the levels of: c_sex 
## Degrees-of-freedom method: kenward-roger 
## P value adjustment: tukey method for comparing a family of 3 estimates

# or for simple comparisons
# pairs(emm, simple = "each")

report(anova(model2b, type = "III"))

## Type 3 ANOVAs only give sensible and informative results when covariates
##   are mean-centered and factors are coded with orthogonal contrasts (such
##   as those produced by `contr.sum`, `contr.poly`, or `contr.helmert`, but
##   *not* by the default `contr.treatment`).

## The ANOVA suggests that:
## 
##   - The main effect of c_age_y_centered is statistically significant and very
## small (F(1) = 4.00, p = 0.046; Eta2 (partial) = 8.77e-03, 95% CI [7.86e-05,
## 1.00])
##   - The main effect of c_sex is statistically not significant and very small
## (F(1) = 0.03, p = 0.857; Eta2 (partial) = 7.17e-05, 95% CI [0.00, 1.00])
##   - The main effect of cbq_fear is statistically not significant and very small
## (F(1) = 1.97, p = 0.161; Eta2 (partial) = 6.95e-03, 95% CI [0.00, 1.00])
##   - The main effect of block is statistically significant and small (F(2) = 3.45,
## p = 0.034; Eta2 (partial) = 0.04, 95% CI [1.64e-03, 1.00])
##   - The main effect of cbq_effortful_control is statistically not significant and
## very small (F(1) = 0.01, p = 0.917; Eta2 (partial) = 3.75e-05, 95% CI [0.00,
## 1.00])
##   - The interaction between cbq_fear and block is statistically not significant
## and small (F(2) = 2.86, p = 0.060; Eta2 (partial) = 0.03, 95% CI [0.00, 1.00])
##   - The interaction between block and cbq_effortful_control is statistically not
## significant and small (F(2) = 2.79, p = 0.064; Eta2 (partial) = 0.03, 95% CI
## [0.00, 1.00])
## 
## Effect sizes were labelled following Field's (2013) recommendations.

Simple intercepts & slopes: EC*Block

interaction_eff2 <- ggpredict(model2, c("cbq_effortful_control", "block"))
interaction_eff2

ggplot(interaction_eff2, aes(x, predicted, color = group)) +
  geom_line(linewidth = 1.2) +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high, fill = group), alpha = 0.2, color = NA) +
  theme_classic() +
  scale_color_manual(
    name = "Block Type",
    values = natparks.pals("Saguaro", 3),
    labels = c(
      "Block 1" = "Baseline",
      "Block 2" = "Stressful",
      "Block 3" = "Recovery"
    )
  ) +
  scale_fill_manual(
    name = "Block Type",
    values = natparks.pals("Saguaro", 3),
    labels = c(
      "Block 1" = "Baseline",
      "Block 2" = "Stressful",
      "Block 3" = "Recovery"
    )
  ) +
  scale_x_continuous(
    name = "Effortful Control (CBQ scores)"
  ) +
  theme(
    axis.text = element_text(size = 19),
    axis.title = element_text(size = 20),
    legend.text = element_text(size = 19),
    legend.title = element_text(face = "bold", size = 19),
    legend.title.align = 0.5
  ) +
  scale_y_continuous(name = "DLPFC Neural Synchrony (predicted)", limits = c(0.29, .42), breaks = seq(0.29, .42, by = 0.02)) # Setting Y-axis limits and ticks

ggsave("figures/aim_2_EC-Block_DLPFC.png", width = 10, height = 8)

Simple intercepts & slopes: F*Block

interaction_eff3 <- ggpredict(model2, c("cbq_fear", "block"))
interaction_eff3

ggplot(interaction_eff3, aes(x, predicted, color = group)) +
  geom_line(linewidth = 1.2) +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high, fill = group), alpha = 0.2, color = NA) +
  theme_classic() +
  scale_color_manual(
    name = "Block Type",
    values = natparks.pals("Saguaro", 3),
    labels = c(
      "Block 1" = "Baseline",
      "Block 2" = "Stressful",
      "Block 3" = "Recovery"
    )
  ) +
  scale_fill_manual(
    name = "Block Type",
    values = natparks.pals("Saguaro", 3),
    labels = c(
      "Block 1" = "Baseline",
      "Block 2" = "Stressful",
      "Block 3" = "Recovery"
    )
  ) +
  scale_x_continuous(
    name = "Fear (CBQ scores)"
  ) +
  theme(
    axis.text = element_text(size = 19),
    axis.title = element_text(size = 20),
    legend.text = element_text(size = 19),
    legend.title = element_text(face = "bold", size = 19),
    legend.title.align = 0.5
  ) +
  scale_y_continuous(name = "DLPFC Neural Synchrony (predicted)", limits = c(0.29, .42), breaks = seq(0.29, .42, by = 0.02)) # Setting Y-axis limits and ticks

ggsave("figures/aim_2_F-Block_DLPFC.png", width = 10, height = 8)

VLPFC

Model

nirs_per_trial_plus_qst$block <- factor(nirs_per_trial_plus_qst$block, levels = c("Block 2", "Block 1", "Block 3"))

model3 <- lmer(
  VLPFC_avg_score ~ c_age_y_centered +
    c_sex + cbq_fear * block + cbq_effortful_control * block + (1 + block | record_id) + (1 | trial),
  data = nirs_per_trial_plus_qst,
  control = lmerControl(
    optimizer = "bobyqa",
    optCtrl = list(maxfun = 5e5)
  )
)

## boundary (singular) fit: see help('isSingular')

summary(model3)

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: VLPFC_avg_score ~ c_age_y_centered + c_sex + cbq_fear * block +  
##     cbq_effortful_control * block + (1 + block | record_id) +      (1 | trial)
##    Data: nirs_per_trial_plus_qst
## Control: lmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 5e+05))
## 
## REML criterion at convergence: -4675.6
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.0718 -0.6666 -0.0149  0.6554  4.0961 
## 
## Random effects:
##  Groups    Name         Variance  Std.Dev. Corr       
##  record_id (Intercept)  7.453e-05 0.008633            
##            blockBlock 1 8.884e-05 0.009426 -0.99      
##            blockBlock 3 2.317e-05 0.004814 -0.90  0.95
##  trial     (Intercept)  0.000e+00 0.000000            
##  Residual               1.451e-03 0.038087            
## Number of obs: 1300, groups:  record_id, 109; trial, 4
## 
## Fixed effects:
##                                      Estimate Std. Error         df t value
## (Intercept)                         3.502e-01  1.826e-02  1.172e+02  19.172
## c_age_y_centered                   -7.134e-04  9.790e-04  1.295e+02  -0.729
## c_sexF                              2.201e-04  2.390e-03  1.301e+02   0.092
## cbq_fear                           -4.273e-03  1.907e-03  1.148e+02  -2.240
## blockBlock 1                        1.599e-02  2.474e-02  1.231e+02   0.646
## blockBlock 3                       -4.570e-03  2.371e-02  1.948e+02  -0.193
## cbq_effortful_control               5.652e-03  3.229e-03  1.195e+02   1.750
## cbq_fear:blockBlock 1               2.637e-03  2.599e-03  1.224e+02   1.015
## cbq_fear:blockBlock 3               5.434e-03  2.503e-03  1.965e+02   2.171
## blockBlock 1:cbq_effortful_control -9.300e-03  4.344e-03  1.232e+02  -2.141
## blockBlock 3:cbq_effortful_control -7.321e-03  4.164e-03  1.948e+02  -1.758
##                                    Pr(>|t|)    
## (Intercept)                          <2e-16 ***
## c_age_y_centered                     0.4675    
## c_sexF                               0.9268    
## cbq_fear                             0.0270 *  
## blockBlock 1                         0.5194    
## blockBlock 3                         0.8474    
## cbq_effortful_control                0.0827 .  
## cbq_fear:blockBlock 1                0.3123    
## cbq_fear:blockBlock 3                0.0311 *  
## blockBlock 1:cbq_effortful_control   0.0343 *  
## blockBlock 3:cbq_effortful_control   0.0803 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) c_g_y_ c_sexF cbq_fr blckB1 blckB3 cbq_f_ cb_:B1 cb_:B3
## c_g_y_cntrd -0.096                                                        
## c_sexF       0.086 -0.251                                                 
## cbq_fear    -0.404  0.048 -0.023                                          
## blockBlock1 -0.739 -0.002 -0.002  0.299                                   
## blockBlock3 -0.696 -0.002 -0.001  0.281  0.521                            
## cbq_ffrtfl_ -0.902  0.101 -0.155 -0.008  0.662  0.623                     
## cbq_fr:blB1  0.297  0.001  0.000 -0.742 -0.402 -0.209  0.011              
## cbq_fr:blB3  0.279  0.002  0.004 -0.696 -0.209 -0.402  0.010  0.518       
## blckB1:cb__  0.666  0.001  0.002  0.011 -0.903 -0.470 -0.733 -0.016 -0.008
## blckB3:cb__  0.628  0.001 -0.001  0.010 -0.470 -0.901 -0.690 -0.008 -0.019
##             bB1:__
## c_g_y_cntrd       
## c_sexF            
## cbq_fear          
## blockBlock1       
## blockBlock3       
## cbq_ffrtfl_       
## cbq_fr:blB1       
## cbq_fr:blB3       
## blckB1:cb__       
## blckB3:cb__  0.520
## optimizer (bobyqa) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')

emm <- emmeans(model3, ~ cbq_effortful_control * block)
pairs(emm)

##  contrast                                                                                     
##  cbq_effortful_control5.17756153846154 Block 2 - cbq_effortful_control5.17756153846154 Block 1
##  cbq_effortful_control5.17756153846154 Block 2 - cbq_effortful_control5.17756153846154 Block 3
##  cbq_effortful_control5.17756153846154 Block 1 - cbq_effortful_control5.17756153846154 Block 3
##   estimate      SE  df t.ratio p.value
##   0.021712 0.00274 106   7.922  <.0001
##   0.020928 0.00263 106   7.957  <.0001
##  -0.000784 0.00263 106  -0.298  0.9523
## 
## Results are averaged over the levels of: c_sex 
## Degrees-of-freedom method: kenward-roger 
## P value adjustment: tukey method for comparing a family of 3 estimates

# or for simple comparisons
# pairs(emm, simple = "each")

emm <- emmeans(model3, ~ cbq_fear * block)
pairs(emm)

##  contrast                                                             estimate
##  cbq_fear3.96526153846154 Block 2 - cbq_fear3.96526153846154 Block 1  0.021712
##  cbq_fear3.96526153846154 Block 2 - cbq_fear3.96526153846154 Block 3  0.020928
##  cbq_fear3.96526153846154 Block 1 - cbq_fear3.96526153846154 Block 3 -0.000784
##       SE  df t.ratio p.value
##  0.00274 106   7.922  <.0001
##  0.00263 106   7.957  <.0001
##  0.00263 106  -0.298  0.9523
## 
## Results are averaged over the levels of: c_sex 
## Degrees-of-freedom method: kenward-roger 
## P value adjustment: tukey method for comparing a family of 3 estimates

report(model3)

## boundary (singular) fit: see help('isSingular')

## Random effect variances not available. Returned R2 does not account for random effects.

## boundary (singular) fit: see help('isSingular')

## Random effect variances not available. Returned R2 does not account for random effects.

## We fitted a linear mixed model (estimated using REML and BOBYQA optimizer) to
## predict VLPFC_avg_score with c_age_y_centered, c_sex, cbq_fear, block and
## cbq_effortful_control (formula: VLPFC_avg_score ~ c_age_y_centered + c_sex +
## cbq_fear * block + cbq_effortful_control * block). The model included block as
## random effects (formula: list(~1 + block | record_id, ~1 | trial)). The model's
## explanatory power related to the fixed effects alone (marginal R2) is 0.07. The
## model's intercept, corresponding to c_age_y_centered = 0, c_sex = M, cbq_fear =
## 0, block = Block 2 and cbq_effortful_control = 0, is at 0.35 (95% CI [0.31,
## 0.39], t(1281) = 19.17, p < .001). Within this model:
## 
##   - The effect of c age y centered is statistically non-significant and negative
## (beta = -7.13e-04, 95% CI [-2.63e-03, 1.21e-03], t(1281) = -0.73, p = 0.466;
## Std. beta = -0.02, 95% CI [-0.08, 0.04])
##   - The effect of c sex [F] is statistically non-significant and positive (beta =
## 2.20e-04, 95% CI [-4.47e-03, 4.91e-03], t(1281) = 0.09, p = 0.927; Std. beta =
## 5.53e-03, 95% CI [-0.11, 0.12])
##   - The effect of cbq fear is statistically significant and negative (beta =
## -4.27e-03, 95% CI [-8.01e-03, -5.31e-04], t(1281) = -2.24, p = 0.025; Std. beta
## = -0.11, 95% CI [-0.21, -0.01])
##   - The effect of block [Block 1] is statistically non-significant and positive
## (beta = 0.02, 95% CI [-0.03, 0.06], t(1281) = 0.65, p = 0.518; Std. beta =
## -0.54, 95% CI [-0.68, -0.41])
##   - The effect of block [Block 3] is statistically non-significant and negative
## (beta = -4.57e-03, 95% CI [-0.05, 0.04], t(1281) = -0.19, p = 0.847; Std. beta
## = -0.53, 95% CI [-0.65, -0.40])
##   - The effect of cbq effortful control is statistically non-significant and
## positive (beta = 5.65e-03, 95% CI [-6.84e-04, 0.01], t(1281) = 1.75, p = 0.080;
## Std. beta = 0.09, 95% CI [-0.01, 0.19])
##   - The effect of cbq fear × block [Block 1] is statistically non-significant and
## positive (beta = 2.64e-03, 95% CI [-2.46e-03, 7.74e-03], t(1281) = 1.01, p =
## 0.311; Std. beta = 0.07, 95% CI [-0.07, 0.20])
##   - The effect of cbq fear × block [Block 3] is statistically significant and
## positive (beta = 5.43e-03, 95% CI [5.24e-04, 0.01], t(1281) = 2.17, p = 0.030;
## Std. beta = 0.14, 95% CI [0.01, 0.27])
##   - The effect of block [Block 1] × cbq effortful control is statistically
## significant and negative (beta = -9.30e-03, 95% CI [-0.02, -7.77e-04], t(1281)
## = -2.14, p = 0.032; Std. beta = -0.15, 95% CI [-0.28, -0.01])
##   - The effect of block [Block 3] × cbq effortful control is statistically
## non-significant and negative (beta = -7.32e-03, 95% CI [-0.02, 8.48e-04],
## t(1281) = -1.76, p = 0.079; Std. beta = -0.12, 95% CI [-0.25, 0.01])
## 
## Standardized parameters were obtained by fitting the model on a standardized
## version of the dataset. 95% Confidence Intervals (CIs) and p-values were
## computed using a Wald t-distribution approximation.

Simple intercepts & slops: EC*Block

interaction_eff4 <- ggpredict(model3, c("cbq_effortful_control", "block"))
interaction_eff4

ggplot(interaction_eff4, aes(x, predicted, color = group)) +
  geom_line(linewidth = 1.2) +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high, fill = group), alpha = 0.2, color = NA) +
  theme_classic() +
  scale_color_manual(
    name = "Block Type",
    values = natparks.pals("Saguaro", 3),
    labels = c(
      "Block 1" = "Baseline",
      "Block 2" = "Stressful",
      "Block 3" = "Recovery"
    )
  ) +
  scale_fill_manual(
    name = "Block Type",
    values = natparks.pals("Saguaro", 3),
    labels = c(
      "Block 1" = "Baseline",
      "Block 2" = "Stressful",
      "Block 3" = "Recovery"
    )
  ) +
  scale_x_continuous(
    name = "Effortful Control (CBQ scores)"
  ) +
  theme(
    axis.text = element_text(size = 19),
    axis.title = element_text(size = 20),
    legend.text = element_text(size = 19),
    legend.title = element_text(face = "bold", size = 19),
    legend.title.align = 0.5
  ) +
  scale_y_continuous(name = "VLPFC Neural Synchrony (predicted)", limits = c(0.29, .42), breaks = seq(0.29, .42, by = 0.02)) # Setting Y-axis limits and ticks

ggsave("figures/aim_2_EC-Block_VLPFC.png.png", width = 10, height = 8)

Simple intercepts & slopes: F*Block

interaction_eff5 <- ggpredict(model3, c("cbq_fear", "block"))
interaction_eff5

ggplot(interaction_eff5, aes(x, predicted, color = group)) +
  geom_line(linewidth = 1.2) +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high, fill = group), alpha = 0.2, color = NA) +
  theme_classic() +
  scale_color_manual(
    name = "Block Type",
    values = natparks.pals("Saguaro", 3),
    labels = c(
      "Block 1" = "Baseline",
      "Block 2" = "Stressful",
      "Block 3" = "Recovery"
    )
  ) +
  scale_fill_manual(
    name = "Block Type",
    values = natparks.pals("Saguaro", 3),
    labels = c(
      "Block 1" = "Baseline",
      "Block 2" = "Stressful",
      "Block 3" = "Recovery"
    )
  ) +
  scale_x_continuous(
    name = "Fear (CBQ scores)"
  ) +
  theme(
    axis.text = element_text(size = 19),
    axis.title = element_text(size = 20),
    legend.text = element_text(size = 19),
    legend.title = element_text(face = "bold", size = 19),
    legend.title.align = 0.5
  ) +
  scale_y_continuous(name = "VLPFC Neural Synchrony (predicted)", limits = c(0.29, .42), breaks = seq(0.29, .42, by = 0.02)) # Setting Y-axis limits and ticks

ggsave("figures/aim_2_F-Block_VLPFC.png", width = 10, height = 8)

Report

tab_model(model1, model2, model3, p.val = "kr", show.df = TRUE, collapse.ci = TRUE, p.style = "stars", digits = 3)

	avg score		DLPFC avg score		VLPFC avg score
Predictors	Estimates	df	Estimates	df	Estimates	df
(Intercept)	0.335 *** (0.298 – 0.371)	106.760	0.314 *** (0.258 – 0.370)	103.316	0.350 *** (0.314 – 0.386)	106.469
c age y centered	0.000 (-0.001 – 0.002)	103.806	0.003 (-0.000 – 0.005)	105.035	-0.001 (-0.003 – 0.001)	103.673
c sex [F]	0.000 (-0.004 – 0.005)	103.822	-0.001 (-0.007 – 0.006)	105.152	0.000 (-0.005 – 0.005)	103.767
cbq fear	-0.000 (-0.004 – 0.004)	105.477	0.007 * (0.001 – 0.013)	103.200	-0.004 * (-0.008 – -0.000)	104.376
block [Block 1]	0.019 (-0.030 – 0.069)	105.863	0.007 (-0.066 – 0.081)	102.631	0.016 (-0.033 – 0.065)	105.716
block [Block 3]	0.029 (-0.015 – 0.074)	105.506	0.084 * (0.008 – 0.159)	103.328	-0.005 (-0.052 – 0.042)	105.308
cbq effortful control	0.006 (-0.001 – 0.012)	108.569	0.005 (-0.005 – 0.015)	105.016	0.006 (-0.001 – 0.012)	108.507
cbq fear × block [Block 1]	-0.000 (-0.005 – 0.005)	105.690	-0.005 (-0.013 – 0.003)	103.303	0.003 (-0.003 – 0.008)	105.065
cbq fear × block [Block 3]	-0.000 (-0.005 – 0.005)	106.514	-0.010 * (-0.018 – -0.002)	105.811	0.005 * (0.000 – 0.010)	105.953
block [Block 1] × cbq effortful control	-0.008 (-0.016 – 0.001)	105.960	-0.002 (-0.014 – 0.011)	102.340	-0.009 * (-0.018 – -0.001)	105.724
block [Block 3] × cbq effortful control	-0.010 * (-0.018 – -0.002)	105.597	-0.013 * (-0.027 – -0.000)	104.111	-0.007 (-0.016 – 0.001)	105.282
Random Effects
σ2	0.00		0.00		0.00
τ00	0.00 record_id		0.00 record_id		0.00 record_id
	0.00 trial		0.00 trial		0.00 trial
τ11	0.00 record_id.blockBlock 1		0.00 record_id.blockBlock 1		0.00 record_id.blockBlock 1
	0.00 record_id.blockBlock 3		0.00 record_id.blockBlock 3		0.00 record_id.blockBlock 3
ρ01	-1.00 record_id.blockBlock 1		-1.00 record_id.blockBlock 1		-0.99 record_id.blockBlock 1
	-1.00 record_id.blockBlock 3		-1.00 record_id.blockBlock 3		-0.90 record_id.blockBlock 3
N	109 record_id		109 record_id		109 record_id
	4 trial		4 trial		4 trial
Observations	1303		1254		1300
Marginal R2 / Conditional R2	0.089 / NA		0.055 / NA		0.074 / NA
p<0.05   ** p<0.01   *** p<0.001

3. To examine the relationship between child temperament, anxiety, and dyadic neural synchrony

Descriptives

nirs_per_trial_plus_qst$block <- factor(nirs_per_trial_plus_qst$block, levels = c("Block 1", "Block 2", "Block 3"))

descriptives_a <- nirs_per_trial_plus_qst %>%
  select(record_id, block, trial, avg_score, DLPFC_avg_score, VLPFC_avg_score) %>%
  group_by(block) %>%
  summarise(
    PFC_mean = round(mean(avg_score, na.rm = TRUE), 2),
    DLPFC_mean = round(mean(DLPFC_avg_score, na.rm = TRUE), 2),
    VLPFC_mean = round(mean(VLPFC_avg_score, na.rm = TRUE), 2),
    PFC_sd = round(sd(avg_score, na.rm = TRUE), 2),
    DLPFC_sd = round(sd(DLPFC_avg_score, na.rm = TRUE), 2),
    VLPFC_sd = round(sd(VLPFC_avg_score, na.rm = TRUE), 2),
    PFC_min = round(min(avg_score, na.rm = TRUE), 2),
    DLPFC_min = round(min(DLPFC_avg_score, na.rm = TRUE), 2),
    VLPFC_min = round(min(VLPFC_avg_score, na.rm = TRUE), 2),
    PFC_max = round(max(avg_score, na.rm = TRUE), 2),
    DLPFC_max = round(max(DLPFC_avg_score, na.rm = TRUE), 2),
    VLPFC_max = round(max(VLPFC_avg_score, na.rm = TRUE), 2)
  ) %>%
  ungroup() %>%
  pivot_longer(PFC_mean:VLPFC_max, names_to = c("ROI", "stat"), names_sep = "_", values_to = "value") %>%
  arrange(ROI) %>%
  unite("name", c("ROI", "block"))

block1 %>%
  select(stai_trait, stai_state, cbq_effortful_control, cbq_fear, anxiety_c) %>%
  pivot_longer(stai_trait:anxiety_c, names_to = "name", values_to = "value") %>%
  group_by(name) %>%
  summarize(
    mean = round(mean(value, na.rm = TRUE), 2),
    sd = round(sd(value, na.rm = TRUE), 2),
    min = round(min(value, na.rm = TRUE), 2),
    max = round(max(value, na.rm = TRUE), 2)
  ) %>%
  pivot_longer(mean:max, names_to = "stat", values_to = "value") %>%
  rbind(descriptives_a) %>%
  pivot_wider(names_from = stat, values_from = value) %>%
  mutate(
    "Mean (SD)" = paste0(mean, " (", sd, ")"),
    Range = paste0(min, "-", max),
    name = as.factor(name)
  ) %>%
  select(name, "Mean (SD)", "Range") %>%
  gt()

name	Mean (SD)	Range
anxiety_c	0.43 (0.26)	0-1.1
cbq_effortful_control	5.18 (0.63)	3.27-6.35
cbq_fear	3.96 (1.06)	1.67-7
stai_state	32.73 (9.17)	18-62
stai_trait	36.62 (8.89)	20-62
DLPFC_Block 1	0.34 (0.05)	0.18-0.49
DLPFC_Block 2	0.37 (0.06)	0.16-0.58
DLPFC_Block 3	0.34 (0.05)	0.2-0.53
PFC_Block 1	0.34 (0.03)	0.24-0.44
PFC_Block 2	0.36 (0.04)	0.27-0.48
PFC_Block 3	0.34 (0.03)	0.27-0.47
VLPFC_Block 1	0.34 (0.04)	0.24-0.48
VLPFC_Block 2	0.36 (0.04)	0.25-0.52
VLPFC_Block 3	0.34 (0.04)	0.25-0.47

Correlations

Table

corr_df <- block1 %>%
  select(stai_trait, stai_state, cbq_effortful_control, cbq_fear, anxiety_c, c_age_m)

# ++++++++++++++++++++++++++++
# flattenCorrMatrix
# ++++++++++++++++++++++++++++
# cormat : matrix of the correlation coefficients
# pmat : matrix of the correlation p-values
flattenCorrMatrix <- function(cormat, pmat) {
  ut <- upper.tri(cormat)
  data.frame(
    row = rownames(cormat)[row(cormat)[ut]],
    column = rownames(cormat)[col(cormat)[ut]],
    cor = round((cormat)[ut], 2),
    p = round(pmat[ut], 6)
  )
}
res2 <- rcorr(as.matrix(corr_df))
res3 <- flattenCorrMatrix(res2$r, res2$P)

res3 %>%
  mutate(test = paste0("r = ", cor, ", p = ", p)) %>%
  select(row, column, test) %>%
  pivot_wider(names_from = column, values_from = test)

Individual

report(cor.test(block1$stai_state, block1$anxiety_c, method = "pearson", conf.level = 0.95))

## Effect sizes were labelled following Funder's (2019) recommendations.
## 
## The Pearson's product-moment correlation between block1$stai_state and
## block1$anxiety_c is positive, statistically significant, and medium (r = 0.29,
## 95% CI [0.20, 0.37], t(432) = 6.22, p < .001)

report(cor.test(block1$stai_state, block1$cbq_effortful_control, method = "pearson", conf.level = 0.95))

## Effect sizes were labelled following Funder's (2019) recommendations.
## 
## The Pearson's product-moment correlation between block1$stai_state and
## block1$cbq_effortful_control is negative, statistically significant, and small
## (r = -0.16, 95% CI [-0.25, -0.07], t(432) = -3.43, p < .001)

report(cor.test(block1$stai_state, block1$cbq_fear, method = "pearson", conf.level = 0.95))

## Effect sizes were labelled following Funder's (2019) recommendations.
## 
## The Pearson's product-moment correlation between block1$stai_state and
## block1$cbq_fear is positive, statistically significant, and small (r = 0.18,
## 95% CI [0.08, 0.27], t(432) = 3.74, p < .001)

report(cor.test(block1$stai_trait, block1$anxiety_c, method = "pearson", conf.level = 0.95))

## Effect sizes were labelled following Funder's (2019) recommendations.
## 
## The Pearson's product-moment correlation between block1$stai_trait and
## block1$anxiety_c is positive, statistically significant, and large (r = 0.32,
## 95% CI [0.24, 0.41], t(432) = 7.10, p < .001)

report(cor.test(block1$stai_trait, block1$cbq_effortful_control, method = "pearson", conf.level = 0.95))

## Effect sizes were labelled following Funder's (2019) recommendations.
## 
## The Pearson's product-moment correlation between block1$stai_trait and
## block1$cbq_effortful_control is negative, statistically not significant, and
## tiny (r = -0.03, 95% CI [-0.12, 0.06], t(432) = -0.62, p = 0.536)

report(cor.test(block1$stai_trait, block1$cbq_fear, method = "pearson", conf.level = 0.95))

## Effect sizes were labelled following Funder's (2019) recommendations.
## 
## The Pearson's product-moment correlation between block1$stai_trait and
## block1$cbq_fear is positive, statistically significant, and very small (r =
## 0.09, 95% CI [1.01e-04, 0.19], t(432) = 1.97, p = 0.050)

report(cor.test(block1$cbq_fear, block1$cbq_effortful_control, method = "pearson", conf.level = 0.95))

## Effect sizes were labelled following Funder's (2019) recommendations.
## 
## The Pearson's product-moment correlation between block1$cbq_fear and
## block1$cbq_effortful_control is positive, statistically not significant, and
## tiny (r = 0.02, 95% CI [-0.08, 0.11], t(432) = 0.36, p = 0.721)

report(cor.test(block1$stai_state, block1$stai_trait, method = "pearson", conf.level = 0.95))

## Effect sizes were labelled following Funder's (2019) recommendations.
## 
## The Pearson's product-moment correlation between block1$stai_state and
## block1$stai_trait is positive, statistically significant, and very large (r =
## 0.78, 95% CI [0.75, 0.82], t(432) = 26.31, p < .001)

report(cor.test(block1$anxiety_c, block1$cbq_fear, method = "pearson", conf.level = 0.95))

## Effect sizes were labelled following Funder's (2019) recommendations.
## 
## The Pearson's product-moment correlation between block1$anxiety_c and
## block1$cbq_fear is positive, statistically significant, and very large (r =
## 0.46, 95% CI [0.39, 0.54], t(432) = 10.90, p < .001)

report(cor.test(block1$anxiety_c, block1$cbq_effortful_control, method = "pearson", conf.level = 0.95))

## Effect sizes were labelled following Funder's (2019) recommendations.
## 
## The Pearson's product-moment correlation between block1$anxiety_c and
## block1$cbq_effortful_control is positive, statistically not significant, and
## tiny (r = 0.04, 95% CI [-0.05, 0.13], t(432) = 0.82, p = 0.413)

PFC

State - Model -block 2 reference

nirs_per_trial_plus_qst$block <- factor(nirs_per_trial_plus_qst$block, levels = c("Block 2", "Block 1", "Block 3"))

model3_1 <- lmer(
  avg_score ~ c_age_y_centered +
    c_sex + cbq_fear * block + cbq_effortful_control * block + anxiety_c * block + stai_state * block + (1 + block | record_id) + (1 | trial),
  data = nirs_per_trial_plus_qst,
  control = lmerControl(
    optimizer = "bobyqa",
    optCtrl = list(maxfun = 5e5)
  )
)

## boundary (singular) fit: see help('isSingular')

summary(model3_1)

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: 
## avg_score ~ c_age_y_centered + c_sex + cbq_fear * block + cbq_effortful_control *  
##     block + anxiety_c * block + stai_state * block + (1 + block |  
##     record_id) + (1 | trial)
##    Data: nirs_per_trial_plus_qst
## Control: lmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 5e+05))
## 
## REML criterion at convergence: -4928
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.9104 -0.6489 -0.0480  0.6275  3.7492 
## 
## Random effects:
##  Groups    Name         Variance  Std.Dev. Corr       
##  record_id (Intercept)  0.0001783 0.01335             
##            blockBlock 1 0.0002928 0.01711  -1.00      
##            blockBlock 3 0.0001259 0.01122  -1.00  1.00
##  trial     (Intercept)  0.0000000 0.00000             
##  Residual               0.0011190 0.03345             
## Number of obs: 1303, groups:  record_id, 109; trial, 4
## 
## Fixed effects:
##                                      Estimate Std. Error         df t value
## (Intercept)                         3.317e-01  2.095e-02  1.094e+02  15.834
## c_age_y_centered                    4.911e-04  8.499e-04  3.217e+02   0.578
## c_sexF                              2.854e-04  2.083e-03  3.232e+02   0.137
## cbq_fear                           -1.090e-04  2.204e-03  1.077e+02  -0.049
## blockBlock 1                        2.357e-02  2.843e-02  1.126e+02   0.829
## blockBlock 3                        2.415e-02  2.549e-02  1.525e+02   0.947
## cbq_effortful_control               5.967e-03  3.346e-03  1.113e+02   1.783
## anxiety_c                          -1.718e-03  9.041e-03  1.075e+02  -0.190
## stai_state                          7.044e-05  2.387e-04  1.078e+02   0.295
## cbq_fear:blockBlock 1              -6.833e-04  3.004e-03  1.125e+02  -0.227
## cbq_fear:blockBlock 3               5.726e-05  2.706e-03  1.554e+02   0.021
## blockBlock 1:cbq_effortful_control -7.946e-03  4.517e-03  1.127e+02  -1.759
## blockBlock 3:cbq_effortful_control -9.370e-03  4.051e-03  1.529e+02  -2.313
## blockBlock 1:anxiety_c              4.901e-03  1.234e-02  1.126e+02   0.397
## blockBlock 3:anxiety_c             -3.392e-03  1.106e-02  1.525e+02  -0.307
## blockBlock 1:stai_state            -9.587e-05  3.253e-04  1.124e+02  -0.295
## blockBlock 3:stai_state             1.256e-04  2.925e-04  1.542e+02   0.430
##                                    Pr(>|t|)    
## (Intercept)                          <2e-16 ***
## c_age_y_centered                     0.5638    
## c_sexF                               0.8911    
## cbq_fear                             0.9606    
## blockBlock 1                         0.4087    
## blockBlock 3                         0.3449    
## cbq_effortful_control                0.0773 .  
## anxiety_c                            0.8496    
## stai_state                           0.7685    
## cbq_fear:blockBlock 1                0.8205    
## cbq_fear:blockBlock 3                0.9831    
## blockBlock 1:cbq_effortful_control   0.0813 .  
## blockBlock 3:cbq_effortful_control   0.0221 *  
## blockBlock 1:anxiety_c               0.6919    
## blockBlock 3:anxiety_c               0.7594    
## blockBlock 1:stai_state              0.7687    
## blockBlock 3:stai_state              0.6681    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## 
## Correlation matrix not shown by default, as p = 17 > 12.
## Use print(x, correlation=TRUE)  or
##     vcov(x)        if you need it

## optimizer (bobyqa) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')

# emm = emmeans(model3_1, ~ cbq_effortful_control*block )
# pairs(emm)
# # or for simple comparisons
# #pairs(emm, simple = "each")

anova(model3_1, type = "III")

report(model3_1)

## boundary (singular) fit: see help('isSingular')

## Random effect variances not available. Returned R2 does not account for random effects.

## boundary (singular) fit: see help('isSingular')

## Random effect variances not available. Returned R2 does not account for random effects.

## We fitted a linear mixed model (estimated using REML and BOBYQA optimizer) to
## predict avg_score with c_age_y_centered, c_sex, cbq_fear, block,
## cbq_effortful_control, anxiety_c and stai_state (formula: avg_score ~
## c_age_y_centered + c_sex + cbq_fear * block + cbq_effortful_control * block +
## anxiety_c * block + stai_state * block). The model included block as random
## effects (formula: list(~1 + block | record_id, ~1 | trial)). The model's
## explanatory power related to the fixed effects alone (marginal R2) is 0.09. The
## model's intercept, corresponding to c_age_y_centered = 0, c_sex = M, cbq_fear =
## 0, block = Block 2, cbq_effortful_control = 0, anxiety_c = 0 and stai_state =
## 0, is at 0.33 (95% CI [0.29, 0.37], t(1278) = 15.83, p < .001). Within this
## model:
## 
##   - The effect of c age y centered is statistically non-significant and positive
## (beta = 4.91e-04, 95% CI [-1.18e-03, 2.16e-03], t(1278) = 0.58, p = 0.563; Std.
## beta = 0.02, 95% CI [-0.04, 0.07])
##   - The effect of c sex [F] is statistically non-significant and positive (beta =
## 2.85e-04, 95% CI [-3.80e-03, 4.37e-03], t(1278) = 0.14, p = 0.891; Std. beta =
## 7.98e-03, 95% CI [-0.11, 0.12])
##   - The effect of cbq fear is statistically non-significant and negative (beta =
## -1.09e-04, 95% CI [-4.43e-03, 4.21e-03], t(1278) = -0.05, p = 0.961; Std. beta
## = -3.20e-03, 95% CI [-0.13, 0.12])
##   - The effect of block [Block 1] is statistically non-significant and positive
## (beta = 0.02, 95% CI [-0.03, 0.08], t(1278) = 0.83, p = 0.407; Std. beta =
## -0.60, 95% CI [-0.75, -0.44])
##   - The effect of block [Block 3] is statistically non-significant and positive
## (beta = 0.02, 95% CI [-0.03, 0.07], t(1278) = 0.95, p = 0.344; Std. beta =
## -0.60, 95% CI [-0.74, -0.46])
##   - The effect of cbq effortful control is statistically non-significant and
## positive (beta = 5.97e-03, 95% CI [-5.98e-04, 0.01], t(1278) = 1.78, p = 0.075;
## Std. beta = 0.11, 95% CI [-0.01, 0.22])
##   - The effect of anxiety c is statistically non-significant and negative (beta =
## -1.72e-03, 95% CI [-0.02, 0.02], t(1278) = -0.19, p = 0.849; Std. beta = -0.01,
## 95% CI [-0.14, 0.12])
##   - The effect of stai state is statistically non-significant and positive (beta
## = 7.04e-05, 95% CI [-3.98e-04, 5.39e-04], t(1278) = 0.30, p = 0.768; Std. beta
## = 0.02, 95% CI [-0.10, 0.14])
##   - The effect of cbq fear × block [Block 1] is statistically non-significant and
## negative (beta = -6.83e-04, 95% CI [-6.58e-03, 5.21e-03], t(1278) = -0.23, p =
## 0.820; Std. beta = -0.02, 95% CI [-0.19, 0.15])
##   - The effect of cbq fear × block [Block 3] is statistically non-significant and
## positive (beta = 5.73e-05, 95% CI [-5.25e-03, 5.37e-03], t(1278) = 0.02, p =
## 0.983; Std. beta = 1.68e-03, 95% CI [-0.15, 0.16])
##   - The effect of block [Block 1] × cbq effortful control is statistically
## non-significant and negative (beta = -7.95e-03, 95% CI [-0.02, 9.15e-04],
## t(1278) = -1.76, p = 0.079; Std. beta = -0.14, 95% CI [-0.30, 0.02])
##   - The effect of block [Block 3] × cbq effortful control is statistically
## significant and negative (beta = -9.37e-03, 95% CI [-0.02, -1.42e-03], t(1278)
## = -2.31, p = 0.021; Std. beta = -0.17, 95% CI [-0.31, -0.03])
##   - The effect of block [Block 1] × anxiety c is statistically non-significant
## and positive (beta = 4.90e-03, 95% CI [-0.02, 0.03], t(1278) = 0.40, p = 0.691;
## Std. beta = 0.04, 95% CI [-0.14, 0.22])
##   - The effect of block [Block 3] × anxiety c is statistically non-significant
## and negative (beta = -3.39e-03, 95% CI [-0.03, 0.02], t(1278) = -0.31, p =
## 0.759; Std. beta = -0.03, 95% CI [-0.19, 0.14])
##   - The effect of block [Block 1] × stai state is statistically non-significant
## and negative (beta = -9.59e-05, 95% CI [-7.34e-04, 5.42e-04], t(1278) = -0.29,
## p = 0.768; Std. beta = -0.02, 95% CI [-0.19, 0.14])
##   - The effect of block [Block 3] × stai state is statistically non-significant
## and positive (beta = 1.26e-04, 95% CI [-4.48e-04, 7.00e-04], t(1278) = 0.43, p
## = 0.668; Std. beta = 0.03, 95% CI [-0.11, 0.18])
## 
## Standardized parameters were obtained by fitting the model on a standardized
## version of the dataset. 95% Confidence Intervals (CIs) and p-values were
## computed using a Wald t-distribution approximation.

model3_1_x <- lmer(
  avg_score ~ c_age_y_centered +
    c_sex + cbq_fear * block + cbq_effortful_control * block + anxiety_c + stai_state * block + (1 + block | record_id) + (1 | trial),
  data = nirs_per_trial_plus_qst,
  control = lmerControl(
    optimizer = "bobyqa",
    optCtrl = list(maxfun = 5e5)
  )
)

## boundary (singular) fit: see help('isSingular')

model3_1_y <- lmer(
  avg_score ~ c_age_y_centered +
    c_sex + cbq_fear * block + cbq_effortful_control * block + anxiety_c + stai_state + (1 + block | record_id) + (1 | trial),
  data = nirs_per_trial_plus_qst,
  control = lmerControl(
    optimizer = "bobyqa",
    optCtrl = list(maxfun = 5e5)
  )
)

## boundary (singular) fit: see help('isSingular')

performance::compare_performance(model3_1, model3_1_x, model3_1_y)

## Random effect variances not available. Returned R2 does not account for random effects.
## Random effect variances not available. Returned R2 does not account for random effects.
## Random effect variances not available. Returned R2 does not account for random effects.

### model_3_1 is the best model <-

Report

tab_model(model3_1, p.val = "kr", show.df = TRUE)

	avg score
Predictors	Estimates	CI	p	df
(Intercept)	0.33	0.29 – 0.37	<0.001	104.96
c age y centered	0.00	-0.00 – 0.00	0.572	101.80
c sex [F]	0.00	-0.00 – 0.00	0.893	101.91
cbq fear	-0.00	-0.00 – 0.00	0.961	103.47
block [Block 1]	0.02	-0.03 – 0.08	0.409	103.85
block [Block 3]	0.02	-0.03 – 0.07	0.346	103.54
cbq effortful control	0.01	-0.00 – 0.01	0.078	106.77
anxiety c	-0.00	-0.02 – 0.02	0.850	103.26
stai state	0.00	-0.00 – 0.00	0.769	103.50
cbq fear × block [Block 1]	-0.00	-0.01 – 0.01	0.821	103.78
cbq fear × block [Block 3]	0.00	-0.01 – 0.01	0.983	104.84
block [Block 1] × cbq effortful control	-0.01	-0.02 – 0.00	0.081	103.96
block [Block 3] × cbq effortful control	-0.01	-0.02 – -0.00	0.023	103.73
block [Block 1] × anxiety c	0.00	-0.02 – 0.03	0.692	103.91
block [Block 3] × anxiety c	-0.00	-0.03 – 0.02	0.760	103.54
block [Block 1] × stai state	-0.00	-0.00 – 0.00	0.769	103.70
block [Block 3] × stai state	0.00	-0.00 – 0.00	0.668	104.31
Random Effects
σ2	0.00
τ00 record_id	0.00
τ00 trial	0.00
τ11 record_id.blockBlock 1	0.00
τ11 record_id.blockBlock 3	0.00
ρ01 record_id.blockBlock 1	-1.00
ρ01 record_id.blockBlock 3	-1.00
N record_id	109
N trial	4
Observations	1303
Marginal R2 / Conditional R2	0.090 / NA

Plot for interaction: EC

interaction_eff3_1 <- ggpredict(model3_1, c("cbq_effortful_control", "block"))
interaction_eff3_1

ggplot(interaction_eff3_1, aes(x, predicted, color = group)) +
  geom_line(linewidth = 1.2) +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high, fill = group), alpha = 0.2, color = NA) +
  theme_classic() +
  scale_color_manual(
    name = "Block Type",
    values = natparks.pals("Saguaro", 3),
    labels = c(
      "Block 1" = "Baseline",
      "Block 2" = "Stressful",
      "Block 3" = "Recovery"
    )
  ) +
  scale_fill_manual(
    name = "Block Type",
    values = natparks.pals("Saguaro", 3),
    labels = c(
      "Block 1" = "Baseline",
      "Block 2" = "Stressful",
      "Block 3" = "Recovery"
    )
  ) +
  scale_x_continuous(
    name = "Effortful Control (CBQ scores)"
  ) +
  theme(
    axis.text = element_text(size = 19),
    axis.title = element_text(size = 20),
    legend.text = element_text(size = 19),
    legend.title = element_text(face = "bold", size = 19),
    legend.title.align = 0.5
  ) +
  scale_y_continuous(name = "PFC Neural Synchrony (predicted)", limits = c(0.29, .42), breaks = seq(0.29, .42, by = 0.02)) # Setting Y-axis limits and ticks

ggsave("figures/aim_3_EC-Block_state_PFC.png", width = 10, height = 8)

Trait - Model -block 2 reference

nirs_per_trial_plus_qst$block <- factor(nirs_per_trial_plus_qst$block, levels = c("Block 2", "Block 1", "Block 3"))

model3_2 <- lmer(
  avg_score ~ c_age_y_centered +
    c_sex + cbq_fear * block + cbq_effortful_control * block + anxiety_c * block + stai_trait * block +
    (1 + block | record_id) + (1 | trial),
  data = nirs_per_trial_plus_qst,
  control = lmerControl(
    optimizer = "bobyqa",
    optCtrl = list(maxfun = 5e5)
  )
)

## boundary (singular) fit: see help('isSingular')

summary(model3_2)

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: 
## avg_score ~ c_age_y_centered + c_sex + cbq_fear * block + cbq_effortful_control *  
##     block + anxiety_c * block + stai_trait * block + (1 + block |  
##     record_id) + (1 | trial)
##    Data: nirs_per_trial_plus_qst
## Control: lmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 5e+05))
## 
## REML criterion at convergence: -4930.2
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.9192 -0.6358 -0.0481  0.6335  3.7537 
## 
## Random effects:
##  Groups    Name         Variance  Std.Dev. Corr       
##  record_id (Intercept)  0.0001719 0.01311             
##            blockBlock 1 0.0002851 0.01688  -1.00      
##            blockBlock 3 0.0001271 0.01127  -1.00  1.00
##  trial     (Intercept)  0.0000000 0.00000             
##  Residual               0.0011187 0.03345             
## Number of obs: 1303, groups:  record_id, 109; trial, 4
## 
## Fixed effects:
##                                      Estimate Std. Error         df t value
## (Intercept)                         3.234e-01  2.057e-02  1.090e+02  15.720
## c_age_y_centered                    4.483e-04  8.460e-04  3.345e+02   0.530
## c_sexF                              2.997e-04  2.067e-03  3.358e+02   0.145
## cbq_fear                            1.033e-04  2.190e-03  1.079e+02   0.047
## blockBlock 1                        3.264e-02  2.802e-02  1.127e+02   1.165
## blockBlock 3                        3.064e-02  2.526e-02  1.500e+02   1.213
## cbq_effortful_control               5.961e-03  3.268e-03  1.112e+02   1.824
## anxiety_c                          -4.712e-03  9.182e-03  1.076e+02  -0.513
## stai_trait                          3.037e-04  2.438e-04  1.077e+02   1.246
## cbq_fear:blockBlock 1              -9.399e-04  2.993e-03  1.127e+02  -0.314
## cbq_fear:blockBlock 3               1.007e-04  2.709e-03  1.525e+02   0.037
## blockBlock 1:cbq_effortful_control -7.915e-03  4.425e-03  1.129e+02  -1.789
## blockBlock 3:cbq_effortful_control -9.715e-03  3.989e-03  1.502e+02  -2.436
## blockBlock 1:anxiety_c              8.211e-03  1.256e-02  1.128e+02   0.654
## blockBlock 3:anxiety_c             -1.753e-03  1.132e-02  1.502e+02  -0.155
## blockBlock 1:stai_trait            -3.492e-04  3.332e-04  1.126e+02  -1.048
## blockBlock 3:stai_trait            -4.031e-05  3.007e-04  1.504e+02  -0.134
##                                    Pr(>|t|)    
## (Intercept)                          <2e-16 ***
## c_age_y_centered                     0.5965    
## c_sexF                               0.8848    
## cbq_fear                             0.9625    
## blockBlock 1                         0.2465    
## blockBlock 3                         0.2270    
## cbq_effortful_control                0.0708 .  
## anxiety_c                            0.6089    
## stai_trait                           0.2156    
## cbq_fear:blockBlock 1                0.7541    
## cbq_fear:blockBlock 3                0.9704    
## blockBlock 1:cbq_effortful_control   0.0763 .  
## blockBlock 3:cbq_effortful_control   0.0160 *  
## blockBlock 1:anxiety_c               0.5147    
## blockBlock 3:anxiety_c               0.8772    
## blockBlock 1:stai_trait              0.2970    
## blockBlock 3:stai_trait              0.8935    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## 
## Correlation matrix not shown by default, as p = 17 > 12.
## Use print(x, correlation=TRUE)  or
##     vcov(x)        if you need it

## optimizer (bobyqa) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')

# emm = emmeans(model3_2, ~ cbq_effortful_control*block )
# pairs(emm)
# # or for simple comparisons
# #pairs(emm, simple = "each")

anova(model3_2, type = "III")

report(model3_2)

## boundary (singular) fit: see help('isSingular')

## Random effect variances not available. Returned R2 does not account for random effects.

## boundary (singular) fit: see help('isSingular')

## Random effect variances not available. Returned R2 does not account for random effects.

## We fitted a linear mixed model (estimated using REML and BOBYQA optimizer) to
## predict avg_score with c_age_y_centered, c_sex, cbq_fear, block,
## cbq_effortful_control, anxiety_c and stai_trait (formula: avg_score ~
## c_age_y_centered + c_sex + cbq_fear * block + cbq_effortful_control * block +
## anxiety_c * block + stai_trait * block). The model included block as random
## effects (formula: list(~1 + block | record_id, ~1 | trial)). The model's
## explanatory power related to the fixed effects alone (marginal R2) is 0.09. The
## model's intercept, corresponding to c_age_y_centered = 0, c_sex = M, cbq_fear =
## 0, block = Block 2, cbq_effortful_control = 0, anxiety_c = 0 and stai_trait =
## 0, is at 0.32 (95% CI [0.28, 0.36], t(1278) = 15.72, p < .001). Within this
## model:
## 
##   - The effect of c age y centered is statistically non-significant and positive
## (beta = 4.48e-04, 95% CI [-1.21e-03, 2.11e-03], t(1278) = 0.53, p = 0.596; Std.
## beta = 0.02, 95% CI [-0.04, 0.07])
##   - The effect of c sex [F] is statistically non-significant and positive (beta =
## 3.00e-04, 95% CI [-3.75e-03, 4.35e-03], t(1278) = 0.15, p = 0.885; Std. beta =
## 8.38e-03, 95% CI [-0.10, 0.12])
##   - The effect of cbq fear is statistically non-significant and positive (beta =
## 1.03e-04, 95% CI [-4.19e-03, 4.40e-03], t(1278) = 0.05, p = 0.962; Std. beta =
## 3.04e-03, 95% CI [-0.12, 0.13])
##   - The effect of block [Block 1] is statistically non-significant and positive
## (beta = 0.03, 95% CI [-0.02, 0.09], t(1278) = 1.16, p = 0.244; Std. beta =
## -0.60, 95% CI [-0.75, -0.44])
##   - The effect of block [Block 3] is statistically non-significant and positive
## (beta = 0.03, 95% CI [-0.02, 0.08], t(1278) = 1.21, p = 0.225; Std. beta =
## -0.60, 95% CI [-0.74, -0.46])
##   - The effect of cbq effortful control is statistically non-significant and
## positive (beta = 5.96e-03, 95% CI [-4.50e-04, 0.01], t(1278) = 1.82, p = 0.068;
## Std. beta = 0.11, 95% CI [-7.94e-03, 0.22])
##   - The effect of anxiety c is statistically non-significant and negative (beta =
## -4.71e-03, 95% CI [-0.02, 0.01], t(1278) = -0.51, p = 0.608; Std. beta = -0.03,
## 95% CI [-0.17, 0.10])
##   - The effect of stai trait is statistically non-significant and positive (beta
## = 3.04e-04, 95% CI [-1.75e-04, 7.82e-04], t(1278) = 1.25, p = 0.213; Std. beta
## = 0.08, 95% CI [-0.04, 0.19])
##   - The effect of cbq fear × block [Block 1] is statistically non-significant and
## negative (beta = -9.40e-04, 95% CI [-6.81e-03, 4.93e-03], t(1278) = -0.31, p =
## 0.754; Std. beta = -0.03, 95% CI [-0.20, 0.15])
##   - The effect of cbq fear × block [Block 3] is statistically non-significant and
## positive (beta = 1.01e-04, 95% CI [-5.21e-03, 5.42e-03], t(1278) = 0.04, p =
## 0.970; Std. beta = 2.96e-03, 95% CI [-0.15, 0.16])
##   - The effect of block [Block 1] × cbq effortful control is statistically
## non-significant and negative (beta = -7.92e-03, 95% CI [-0.02, 7.66e-04],
## t(1278) = -1.79, p = 0.074; Std. beta = -0.14, 95% CI [-0.29, 0.01])
##   - The effect of block [Block 3] × cbq effortful control is statistically
## significant and negative (beta = -9.72e-03, 95% CI [-0.02, -1.89e-03], t(1278)
## = -2.44, p = 0.015; Std. beta = -0.17, 95% CI [-0.31, -0.03])
##   - The effect of block [Block 1] × anxiety c is statistically non-significant
## and positive (beta = 8.21e-03, 95% CI [-0.02, 0.03], t(1278) = 0.65, p = 0.513;
## Std. beta = 0.06, 95% CI [-0.12, 0.24])
##   - The effect of block [Block 3] × anxiety c is statistically non-significant
## and negative (beta = -1.75e-03, 95% CI [-0.02, 0.02], t(1278) = -0.15, p =
## 0.877; Std. beta = -0.01, 95% CI [-0.18, 0.15])
##   - The effect of block [Block 1] × stai trait is statistically non-significant
## and negative (beta = -3.49e-04, 95% CI [-1.00e-03, 3.05e-04], t(1278) = -1.05,
## p = 0.295; Std. beta = -0.09, 95% CI [-0.25, 0.08])
##   - The effect of block [Block 3] × stai trait is statistically non-significant
## and negative (beta = -4.03e-05, 95% CI [-6.30e-04, 5.50e-04], t(1278) = -0.13,
## p = 0.893; Std. beta = -1.00e-02, 95% CI [-0.16, 0.14])
## 
## Standardized parameters were obtained by fitting the model on a standardized
## version of the dataset. 95% Confidence Intervals (CIs) and p-values were
## computed using a Wald t-distribution approximation.

performance::check_outliers(model3_2)

## OK: No outliers detected.
## - Based on the following method and threshold: cook (0.907).
## - For variable: (Whole model)

model3_2_x <- lmer(
  avg_score ~ c_age_y_centered +
    c_sex + cbq_fear * block + cbq_effortful_control * block + anxiety_c + stai_trait * block +
    (1 + block | record_id) + (1 | trial),
  data = nirs_per_trial_plus_qst,
  control = lmerControl(
    optimizer = "bobyqa",
    optCtrl = list(maxfun = 5e5)
  )
)

## boundary (singular) fit: see help('isSingular')

model3_2_y <- lmer(
  avg_score ~ c_age_y_centered +
    c_sex + cbq_fear * block + cbq_effortful_control * block + anxiety_c + stai_trait +
    (1 + block | record_id) + (1 | trial),
  data = nirs_per_trial_plus_qst,
  control = lmerControl(
    optimizer = "bobyqa",
    optCtrl = list(maxfun = 5e5)
  )
)

## boundary (singular) fit: see help('isSingular')

performance::compare_performance(model3_2, model3_2_x, model3_2_y)

## Random effect variances not available. Returned R2 does not account for random effects.
## Random effect variances not available. Returned R2 does not account for random effects.
## Random effect variances not available. Returned R2 does not account for random effects.

Plot for interaction: EC

interaction_eff3_2 <- ggpredict(model3_2, c("cbq_effortful_control", "block"))
interaction_eff3_2

ggplot(interaction_eff3_2, aes(x, predicted, color = group)) +
  geom_line(linewidth = 1.2) +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high, fill = group), alpha = 0.2, color = NA) +
  theme_classic() +
  scale_color_manual(
    name = "Block Type",
    values = natparks.pals("Saguaro", 3),
    labels = c(
      "Block 1" = "Baseline",
      "Block 2" = "Stressful",
      "Block 3" = "Recovery"
    )
  ) +
  scale_fill_manual(
    name = "Block Type",
    values = natparks.pals("Saguaro", 3),
    labels = c(
      "Block 1" = "Baseline",
      "Block 2" = "Stressful",
      "Block 3" = "Recovery"
    )
  ) +
  scale_x_continuous(
    name = "Effortful Control (CBQ scores)"
  ) +
  theme(
    axis.text = element_text(size = 19),
    axis.title = element_text(size = 20),
    legend.text = element_text(size = 19),
    legend.title = element_text(face = "bold", size = 19),
    legend.title.align = 0.5
  ) +
  scale_y_continuous(name = "PFC Neural Synchrony (predicted)", limits = c(0.29, .42), breaks = seq(0.29, .42, by = 0.02)) # Setting Y-axis limits and ticks

ggsave("figures/aim_3_EC-Block_trait_PFC.png", width = 10, height = 8)

DLPFC

State - Model -block 2 reference

nirs_per_trial_plus_qst$block <- factor(nirs_per_trial_plus_qst$block, levels = c("Block 2", "Block 1", "Block 3"))

model3_3 <- lmer(
  DLPFC_avg_score ~ c_age_y_centered +
    c_sex + cbq_fear * block + cbq_effortful_control * block + anxiety_c * block + stai_state * block
    + (1 + block | record_id) + (1 | trial),
  data = nirs_per_trial_plus_qst,
  control = lmerControl(
    optimizer = "bobyqa",
    optCtrl = list(maxfun = 5e5)
  )
)

## boundary (singular) fit: see help('isSingular')

summary(model3_3)

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: DLPFC_avg_score ~ c_age_y_centered + c_sex + cbq_fear * block +  
##     cbq_effortful_control * block + anxiety_c * block + stai_state *  
##     block + (1 + block | record_id) + (1 | trial)
##    Data: nirs_per_trial_plus_qst
## Control: lmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 5e+05))
## 
## REML criterion at convergence: -3606.9
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.2594 -0.6981 -0.0131  0.6264  3.6368 
## 
## Random effects:
##  Groups    Name         Variance  Std.Dev. Corr       
##  record_id (Intercept)  0.0003510 0.01874             
##            blockBlock 1 0.0003852 0.01963  -1.00      
##            blockBlock 3 0.0005459 0.02337  -1.00  1.00
##  trial     (Intercept)  0.0000000 0.00000             
##  Residual               0.0028138 0.05305             
## Number of obs: 1254, groups:  record_id, 109; trial, 4
## 
## Fixed effects:
##                                      Estimate Std. Error         df t value
## (Intercept)                         3.062e-01  3.210e-02  1.012e+02   9.538
## c_age_y_centered                    2.736e-03  1.350e-03  4.470e+02   2.027
## c_sexF                             -9.361e-04  3.312e-03  4.478e+02  -0.283
## cbq_fear                            7.544e-03  3.406e-03  1.008e+02   2.215
## blockBlock 1                        1.506e-02  4.161e-02  1.257e+02   0.362
## blockBlock 3                        7.718e-02  4.342e-02  1.094e+02   1.778
## cbq_effortful_control               5.302e-03  5.128e-03  1.031e+02   1.034
## anxiety_c                          -6.767e-03  1.398e-02  1.014e+02  -0.484
## stai_state                          1.770e-04  3.662e-04  9.997e+01   0.483
## cbq_fear:blockBlock 1              -7.124e-03  4.417e-03  1.265e+02  -1.613
## cbq_fear:blockBlock 3              -1.149e-02  4.637e-03  1.129e+02  -2.478
## blockBlock 1:cbq_effortful_control -2.132e-03  6.614e-03  1.263e+02  -0.322
## blockBlock 3:cbq_effortful_control -1.311e-02  6.912e-03  1.104e+02  -1.897
## blockBlock 1:anxiety_c              1.947e-02  1.820e-02  1.288e+02   1.069
## blockBlock 3:anxiety_c              1.241e-02  1.892e-02  1.107e+02   0.656
## blockBlock 1:stai_state            -1.581e-04  4.767e-04  1.258e+02  -0.332
## blockBlock 3:stai_state             1.955e-04  4.978e-04  1.099e+02   0.393
##                                    Pr(>|t|)    
## (Intercept)                        9.24e-16 ***
## c_age_y_centered                     0.0433 *  
## c_sexF                               0.7776    
## cbq_fear                             0.0290 *  
## blockBlock 1                         0.7180    
## blockBlock 3                         0.0782 .  
## cbq_effortful_control                0.3035    
## anxiety_c                            0.6294    
## stai_state                           0.6300    
## cbq_fear:blockBlock 1                0.1093    
## cbq_fear:blockBlock 3                0.0147 *  
## blockBlock 1:cbq_effortful_control   0.7477    
## blockBlock 3:cbq_effortful_control   0.0605 .  
## blockBlock 1:anxiety_c               0.2869    
## blockBlock 3:anxiety_c               0.5135    
## blockBlock 1:stai_state              0.7408    
## blockBlock 3:stai_state              0.6953    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## 
## Correlation matrix not shown by default, as p = 17 > 12.
## Use print(x, correlation=TRUE)  or
##     vcov(x)        if you need it

## optimizer (bobyqa) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')

# emm = emmeans(model3_3, ~ cbq_effortful_control*block )
# pairs(emm)
# # or for simple comparisons
# #pairs(emm, simple = "each")

anova(model3_3, type = "III")

report(model3_3)

## boundary (singular) fit: see help('isSingular')

## Random effect variances not available. Returned R2 does not account for random effects.

## boundary (singular) fit: see help('isSingular')

## Random effect variances not available. Returned R2 does not account for random effects.

## We fitted a linear mixed model (estimated using REML and BOBYQA optimizer) to
## predict DLPFC_avg_score with c_age_y_centered, c_sex, cbq_fear, block,
## cbq_effortful_control, anxiety_c and stai_state (formula: DLPFC_avg_score ~
## c_age_y_centered + c_sex + cbq_fear * block + cbq_effortful_control * block +
## anxiety_c * block + stai_state * block). The model included block as random
## effects (formula: list(~1 + block | record_id, ~1 | trial)). The model's
## explanatory power related to the fixed effects alone (marginal R2) is 0.06. The
## model's intercept, corresponding to c_age_y_centered = 0, c_sex = M, cbq_fear =
## 0, block = Block 2, cbq_effortful_control = 0, anxiety_c = 0 and stai_state =
## 0, is at 0.31 (95% CI [0.24, 0.37], t(1229) = 9.54, p < .001). Within this
## model:
## 
##   - The effect of c age y centered is statistically significant and positive
## (beta = 2.74e-03, 95% CI [8.74e-05, 5.39e-03], t(1229) = 2.03, p = 0.043; Std.
## beta = 0.06, 95% CI [1.90e-03, 0.12])
##   - The effect of c sex [F] is statistically non-significant and negative (beta =
## -9.36e-04, 95% CI [-7.43e-03, 5.56e-03], t(1229) = -0.28, p = 0.777; Std. beta
## = -0.02, 95% CI [-0.13, 0.10])
##   - The effect of cbq fear is statistically significant and positive (beta =
## 7.54e-03, 95% CI [8.61e-04, 0.01], t(1229) = 2.21, p = 0.027; Std. beta = 0.14,
## 95% CI [0.02, 0.27])
##   - The effect of block [Block 1] is statistically non-significant and positive
## (beta = 0.02, 95% CI [-0.07, 0.10], t(1229) = 0.36, p = 0.717; Std. beta =
## -0.38, 95% CI [-0.53, -0.23])
##   - The effect of block [Block 3] is statistically non-significant and positive
## (beta = 0.08, 95% CI [-8.00e-03, 0.16], t(1229) = 1.78, p = 0.076; Std. beta =
## -0.44, 95% CI [-0.59, -0.29])
##   - The effect of cbq effortful control is statistically non-significant and
## positive (beta = 5.30e-03, 95% CI [-4.76e-03, 0.02], t(1229) = 1.03, p = 0.301;
## Std. beta = 0.06, 95% CI [-0.05, 0.18])
##   - The effect of anxiety c is statistically non-significant and negative (beta =
## -6.77e-03, 95% CI [-0.03, 0.02], t(1229) = -0.48, p = 0.628; Std. beta = -0.03,
## 95% CI [-0.16, 0.10])
##   - The effect of stai state is statistically non-significant and positive (beta
## = 1.77e-04, 95% CI [-5.41e-04, 8.95e-04], t(1229) = 0.48, p = 0.629; Std. beta
## = 0.03, 95% CI [-0.09, 0.15])
##   - The effect of cbq fear × block [Block 1] is statistically non-significant and
## negative (beta = -7.12e-03, 95% CI [-0.02, 1.54e-03], t(1229) = -1.61, p =
## 0.107; Std. beta = -0.14, 95% CI [-0.30, 0.03])
##   - The effect of cbq fear × block [Block 3] is statistically significant and
## negative (beta = -0.01, 95% CI [-0.02, -2.39e-03], t(1229) = -2.48, p = 0.013;
## Std. beta = -0.22, 95% CI [-0.39, -0.05])
##   - The effect of block [Block 1] × cbq effortful control is statistically
## non-significant and negative (beta = -2.13e-03, 95% CI [-0.02, 0.01], t(1229) =
## -0.32, p = 0.747; Std. beta = -0.02, 95% CI [-0.17, 0.12])
##   - The effect of block [Block 3] × cbq effortful control is statistically
## non-significant and negative (beta = -0.01, 95% CI [-0.03, 4.51e-04], t(1229) =
## -1.90, p = 0.058; Std. beta = -0.15, 95% CI [-0.31, 5.17e-03])
##   - The effect of block [Block 1] × anxiety c is statistically non-significant
## and positive (beta = 0.02, 95% CI [-0.02, 0.06], t(1229) = 1.07, p = 0.285;
## Std. beta = 0.09, 95% CI [-0.08, 0.26])
##   - The effect of block [Block 3] × anxiety c is statistically non-significant
## and positive (beta = 0.01, 95% CI [-0.02, 0.05], t(1229) = 0.66, p = 0.512;
## Std. beta = 0.06, 95% CI [-0.12, 0.24])
##   - The effect of block [Block 1] × stai state is statistically non-significant
## and negative (beta = -1.58e-04, 95% CI [-1.09e-03, 7.77e-04], t(1229) = -0.33,
## p = 0.740; Std. beta = -0.03, 95% CI [-0.18, 0.13])
##   - The effect of block [Block 3] × stai state is statistically non-significant
## and positive (beta = 1.95e-04, 95% CI [-7.81e-04, 1.17e-03], t(1229) = 0.39, p
## = 0.695; Std. beta = 0.03, 95% CI [-0.13, 0.20])
## 
## Standardized parameters were obtained by fitting the model on a standardized
## version of the dataset. 95% Confidence Intervals (CIs) and p-values were
## computed using a Wald t-distribution approximation.

model3_3_x <- lmer(
  DLPFC_avg_score ~ c_age_y_centered +
    c_sex + cbq_fear * block + cbq_effortful_control * block + anxiety_c + stai_state * block +
    (1 + block | record_id) + (1 | trial),
  data = nirs_per_trial_plus_qst,
  control = lmerControl(
    optimizer = "bobyqa",
    optCtrl = list(maxfun = 5e5)
  )
)

## boundary (singular) fit: see help('isSingular')

model3_3_y <- lmer(
  DLPFC_avg_score ~ c_age_y_centered +
    c_sex + cbq_fear * block + cbq_effortful_control * block + anxiety_c + stai_state +
    (1 + block | record_id) + (1 | trial),
  data = nirs_per_trial_plus_qst,
  control = lmerControl(
    optimizer = "bobyqa",
    optCtrl = list(maxfun = 5e5)
  )
)

## boundary (singular) fit: see help('isSingular')

performance::compare_performance(model3_3, model3_3_x, model3_3_y)

## Random effect variances not available. Returned R2 does not account for random effects.
## Random effect variances not available. Returned R2 does not account for random effects.
## Random effect variances not available. Returned R2 does not account for random effects.

## best model is model3_3

Plot for interaction: Fear

interaction_eff3_3 <- ggpredict(model3_3, c("cbq_fear", "block"))
interaction_eff3_3

ggplot(interaction_eff3_3, aes(x, predicted, color = group)) +
  geom_line(linewidth = 1.2) +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high, fill = group), alpha = 0.2, color = NA) +
  theme_classic() +
  scale_color_manual(
    name = "Block Type",
    values = natparks.pals("Saguaro", 3),
    labels = c(
      "Block 1" = "Baseline",
      "Block 2" = "Stressful",
      "Block 3" = "Recovery"
    )
  ) +
  scale_fill_manual(
    name = "Block Type",
    values = natparks.pals("Saguaro", 3),
    labels = c(
      "Block 1" = "Baseline",
      "Block 2" = "Stressful",
      "Block 3" = "Recovery"
    )
  ) +
  scale_x_continuous(
    name = "Fear (CBQ scores)"
  ) +
  theme(
    axis.text = element_text(size = 19),
    axis.title = element_text(size = 20),
    legend.text = element_text(size = 19),
    legend.title = element_text(face = "bold", size = 19),
    legend.title.align = 0.5
  ) +
  scale_y_continuous(name = "DLPFC Neural Synchrony (predicted)", limits = c(0.29, .42), breaks = seq(0.29, .42, by = 0.02)) # Setting Y-axis limits and ticks

ggsave("figures/aim_3_F-Block_state_DLPFC.png", width = 10, height = 8)

Trait - Model -block 2 reference

nirs_per_trial_plus_qst$block <- factor(nirs_per_trial_plus_qst$block, levels = c("Block 2", "Block 1", "Block 3"))

model3_4 <- lmer(
  DLPFC_avg_score ~ c_age_y_centered +
    c_sex + cbq_fear * block + cbq_effortful_control * block + anxiety_c * block + stai_trait * block + (1 + block | record_id) + (1 | trial),
  data = nirs_per_trial_plus_qst,
  control = lmerControl(
    optimizer = "bobyqa",
    optCtrl = list(maxfun = 5e5)
  )
)

## boundary (singular) fit: see help('isSingular')

summary(model3_4)

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: DLPFC_avg_score ~ c_age_y_centered + c_sex + cbq_fear * block +  
##     cbq_effortful_control * block + anxiety_c * block + stai_trait *  
##     block + (1 + block | record_id) + (1 | trial)
##    Data: nirs_per_trial_plus_qst
## Control: lmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 5e+05))
## 
## REML criterion at convergence: -3608.9
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.2215 -0.6979 -0.0189  0.6367  3.6747 
## 
## Random effects:
##  Groups    Name         Variance  Std.Dev. Corr       
##  record_id (Intercept)  0.0003336 0.01826             
##            blockBlock 1 0.0003708 0.01926  -1.00      
##            blockBlock 3 0.0005398 0.02323  -1.00  1.00
##  trial     (Intercept)  0.0000000 0.00000             
##  Residual               0.0028126 0.05303             
## Number of obs: 1254, groups:  record_id, 109; trial, 4
## 
## Fixed effects:
##                                      Estimate Std. Error         df t value
## (Intercept)                         2.949e-01  3.151e-02  1.015e+02   9.359
## c_age_y_centered                    2.634e-03  1.344e-03  4.777e+02   1.959
## c_sexF                             -8.370e-04  3.287e-03  4.793e+02  -0.255
## cbq_fear                            7.895e-03  3.378e-03  1.012e+02   2.337
## blockBlock 1                        2.655e-02  4.108e-02  1.270e+02   0.646
## blockBlock 3                        8.921e-02  4.297e-02  1.093e+02   2.076
## cbq_effortful_control               5.109e-03  5.002e-03  1.034e+02   1.021
## anxiety_c                          -1.094e-02  1.415e-02  1.015e+02  -0.773
## stai_trait                          5.050e-04  3.750e-04  1.005e+02   1.347
## cbq_fear:blockBlock 1              -7.464e-03  4.399e-03  1.271e+02  -1.697
## cbq_fear:blockBlock 3              -1.139e-02  4.630e-03  1.122e+02  -2.461
## blockBlock 1:cbq_effortful_control -2.014e-03  6.488e-03  1.274e+02  -0.310
## blockBlock 3:cbq_effortful_control -1.366e-02  6.791e-03  1.099e+02  -2.012
## blockBlock 1:anxiety_c              2.358e-02  1.850e-02  1.289e+02   1.274
## blockBlock 3:anxiety_c              1.532e-02  1.932e-02  1.102e+02   0.793
## blockBlock 1:stai_trait            -4.842e-04  4.886e-04  1.263e+02  -0.991
## blockBlock 3:stai_trait            -1.214e-04  5.120e-04  1.095e+02  -0.237
##                                    Pr(>|t|)    
## (Intercept)                        2.22e-15 ***
## c_age_y_centered                     0.0507 .  
## c_sexF                               0.7991    
## cbq_fear                             0.0214 *  
## blockBlock 1                         0.5193    
## blockBlock 3                         0.0402 *  
## cbq_effortful_control                0.3095    
## anxiety_c                            0.4414    
## stai_trait                           0.1811    
## cbq_fear:blockBlock 1                0.0922 .  
## cbq_fear:blockBlock 3                0.0154 *  
## blockBlock 1:cbq_effortful_control   0.7567    
## blockBlock 3:cbq_effortful_control   0.0467 *  
## blockBlock 1:anxiety_c               0.2049    
## blockBlock 3:anxiety_c               0.4296    
## blockBlock 1:stai_trait              0.3236    
## blockBlock 3:stai_trait              0.8131    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## 
## Correlation matrix not shown by default, as p = 17 > 12.
## Use print(x, correlation=TRUE)  or
##     vcov(x)        if you need it

## optimizer (bobyqa) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')

# emm = emmeans(model3_4, ~ cbq_effortful_control*block )
# pairs(emm)
# # or for simple comparisons
# #pairs(emm, simple = "each")

anova(model3_4, type = "III")

report(model3_4)

## boundary (singular) fit: see help('isSingular')

## Random effect variances not available. Returned R2 does not account for random effects.

## boundary (singular) fit: see help('isSingular')

## Random effect variances not available. Returned R2 does not account for random effects.

## We fitted a linear mixed model (estimated using REML and BOBYQA optimizer) to
## predict DLPFC_avg_score with c_age_y_centered, c_sex, cbq_fear, block,
## cbq_effortful_control, anxiety_c and stai_trait (formula: DLPFC_avg_score ~
## c_age_y_centered + c_sex + cbq_fear * block + cbq_effortful_control * block +
## anxiety_c * block + stai_trait * block). The model included block as random
## effects (formula: list(~1 + block | record_id, ~1 | trial)). The model's
## explanatory power related to the fixed effects alone (marginal R2) is 0.06. The
## model's intercept, corresponding to c_age_y_centered = 0, c_sex = M, cbq_fear =
## 0, block = Block 2, cbq_effortful_control = 0, anxiety_c = 0 and stai_trait =
## 0, is at 0.29 (95% CI [0.23, 0.36], t(1229) = 9.36, p < .001). Within this
## model:
## 
##   - The effect of c age y centered is statistically non-significant and positive
## (beta = 2.63e-03, 95% CI [-3.61e-06, 5.27e-03], t(1229) = 1.96, p = 0.050; Std.
## beta = 0.06, 95% CI [-7.83e-05, 0.11])
##   - The effect of c sex [F] is statistically non-significant and negative (beta =
## -8.37e-04, 95% CI [-7.29e-03, 5.61e-03], t(1229) = -0.25, p = 0.799; Std. beta
## = -0.02, 95% CI [-0.13, 0.10])
##   - The effect of cbq fear is statistically significant and positive (beta =
## 7.89e-03, 95% CI [1.27e-03, 0.01], t(1229) = 2.34, p = 0.020; Std. beta = 0.15,
## 95% CI [0.02, 0.28])
##   - The effect of block [Block 1] is statistically non-significant and positive
## (beta = 0.03, 95% CI [-0.05, 0.11], t(1229) = 0.65, p = 0.518; Std. beta =
## -0.38, 95% CI [-0.53, -0.23])
##   - The effect of block [Block 3] is statistically significant and positive (beta
## = 0.09, 95% CI [4.91e-03, 0.17], t(1229) = 2.08, p = 0.038; Std. beta = -0.44,
## 95% CI [-0.59, -0.29])
##   - The effect of cbq effortful control is statistically non-significant and
## positive (beta = 5.11e-03, 95% CI [-4.70e-03, 0.01], t(1229) = 1.02, p = 0.307;
## Std. beta = 0.06, 95% CI [-0.05, 0.17])
##   - The effect of anxiety c is statistically non-significant and negative (beta =
## -0.01, 95% CI [-0.04, 0.02], t(1229) = -0.77, p = 0.440; Std. beta = -0.05, 95%
## CI [-0.18, 0.08])
##   - The effect of stai trait is statistically non-significant and positive (beta
## = 5.05e-04, 95% CI [-2.31e-04, 1.24e-03], t(1229) = 1.35, p = 0.178; Std. beta
## = 0.08, 95% CI [-0.04, 0.20])
##   - The effect of cbq fear × block [Block 1] is statistically non-significant and
## negative (beta = -7.46e-03, 95% CI [-0.02, 1.17e-03], t(1229) = -1.70, p =
## 0.090; Std. beta = -0.14, 95% CI [-0.31, 0.02])
##   - The effect of cbq fear × block [Block 3] is statistically significant and
## negative (beta = -0.01, 95% CI [-0.02, -2.31e-03], t(1229) = -2.46, p = 0.014;
## Std. beta = -0.22, 95% CI [-0.39, -0.04])
##   - The effect of block [Block 1] × cbq effortful control is statistically
## non-significant and negative (beta = -2.01e-03, 95% CI [-0.01, 0.01], t(1229) =
## -0.31, p = 0.756; Std. beta = -0.02, 95% CI [-0.17, 0.12])
##   - The effect of block [Block 3] × cbq effortful control is statistically
## significant and negative (beta = -0.01, 95% CI [-0.03, -3.40e-04], t(1229) =
## -2.01, p = 0.044; Std. beta = -0.16, 95% CI [-0.31, -3.90e-03])
##   - The effect of block [Block 1] × anxiety c is statistically non-significant
## and positive (beta = 0.02, 95% CI [-0.01, 0.06], t(1229) = 1.27, p = 0.203;
## Std. beta = 0.11, 95% CI [-0.06, 0.29])
##   - The effect of block [Block 3] × anxiety c is statistically non-significant
## and positive (beta = 0.02, 95% CI [-0.02, 0.05], t(1229) = 0.79, p = 0.428;
## Std. beta = 0.07, 95% CI [-0.11, 0.25])
##   - The effect of block [Block 1] × stai trait is statistically non-significant
## and negative (beta = -4.84e-04, 95% CI [-1.44e-03, 4.74e-04], t(1229) = -0.99,
## p = 0.322; Std. beta = -0.08, 95% CI [-0.23, 0.08])
##   - The effect of block [Block 3] × stai trait is statistically non-significant
## and negative (beta = -1.21e-04, 95% CI [-1.13e-03, 8.83e-04], t(1229) = -0.24,
## p = 0.813; Std. beta = -0.02, 95% CI [-0.18, 0.14])
## 
## Standardized parameters were obtained by fitting the model on a standardized
## version of the dataset. 95% Confidence Intervals (CIs) and p-values were
## computed using a Wald t-distribution approximation.

performance::check_outliers(model3_4)

## OK: No outliers detected.
## - Based on the following method and threshold: cook (0.907).
## - For variable: (Whole model)

model3_4_x <- lmer(
  DLPFC_avg_score ~ c_age_y_centered +
    c_sex + cbq_fear * block + cbq_effortful_control * block + anxiety_c + stai_trait * block +
    (1 + block | record_id) + (1 | trial),
  data = nirs_per_trial_plus_qst,
  control = lmerControl(
    optimizer = "bobyqa",
    optCtrl = list(maxfun = 5e5)
  )
)

## boundary (singular) fit: see help('isSingular')

model3_4_y <- lmer(
  DLPFC_avg_score ~ c_age_y_centered +
    c_sex + cbq_fear * block + cbq_effortful_control * block + anxiety_c + stai_trait +
    (1 + block | record_id) + (1 | trial),
  data = nirs_per_trial_plus_qst,
  control = lmerControl(
    optimizer = "bobyqa",
    optCtrl = list(maxfun = 5e5)
  )
)

## boundary (singular) fit: see help('isSingular')

performance::compare_performance(model3_4, model3_4_x, model3_4_y)

## Random effect variances not available. Returned R2 does not account for random effects.
## Random effect variances not available. Returned R2 does not account for random effects.
## Random effect variances not available. Returned R2 does not account for random effects.

Trait - Model -block 3 reference

nirs_per_trial_plus_qst$block <- factor(nirs_per_trial_plus_qst$block, levels = c("Block 3", "Block 1", "Block 2"))

model3_4b <- lmer(
  DLPFC_avg_score ~ c_age_y_centered +
    c_sex + cbq_fear * block + cbq_effortful_control * block + anxiety_c * block + stai_trait * block + (1 + block | record_id) + (1 | trial),
  data = nirs_per_trial_plus_qst,
  control = lmerControl(
    optimizer = "bobyqa",
    optCtrl = list(maxfun = 5e5)
  )
)

## boundary (singular) fit: see help('isSingular')

## Warning: Model failed to converge with 1 negative eigenvalue: -1.8e-02

summary(model3_4b)

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: DLPFC_avg_score ~ c_age_y_centered + c_sex + cbq_fear * block +  
##     cbq_effortful_control * block + anxiety_c * block + stai_trait *  
##     block + (1 + block | record_id) + (1 | trial)
##    Data: nirs_per_trial_plus_qst
## Control: lmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 5e+05))
## 
## REML criterion at convergence: -3607.5
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.2508 -0.7004 -0.0217  0.6417  3.6584 
## 
## Random effects:
##  Groups    Name         Variance  Std.Dev. Corr       
##  record_id (Intercept)  0.000e+00 0.000000            
##            blockBlock 1 1.306e-06 0.001143   NaN      
##            blockBlock 2 3.227e-04 0.017964   NaN -1.00
##  trial     (Intercept)  0.000e+00 0.000000            
##  Residual               2.824e-03 0.053143            
## Number of obs: 1254, groups:  record_id, 109; trial, 4
## 
## Fixed effects:
##                                      Estimate Std. Error         df t value
## (Intercept)                         3.835e-01  2.611e-02  1.127e+03  14.686
## c_age_y_centered                    2.656e-03  1.378e-03  2.812e+02   1.927
## c_sexF                             -9.466e-04  3.370e-03  2.822e+02  -0.281
## cbq_fear                           -3.373e-03  2.811e-03  1.124e+03  -1.200
## blockBlock 1                       -6.215e-02  3.667e-02  1.042e+03  -1.695
## blockBlock 2                       -8.866e-02  4.062e-02  1.189e+02  -2.183
## cbq_effortful_control              -8.503e-03  4.181e-03  1.129e+03  -2.034
## anxiety_c                           4.170e-03  1.162e-02  1.123e+03   0.359
## stai_trait                          3.833e-04  3.082e-04  1.124e+03   1.244
## cbq_fear:blockBlock 1               3.806e-03  3.935e-03  1.043e+03   0.967
## cbq_fear:blockBlock 2               1.123e-02  4.381e-03  1.223e+02   2.562
## blockBlock 1:cbq_effortful_control  1.162e-02  5.809e-03  1.043e+03   2.001
## blockBlock 2:cbq_effortful_control  1.366e-02  6.420e-03  1.196e+02   2.128
## blockBlock 1:anxiety_c              8.457e-03  1.646e-02  1.042e+03   0.514
## blockBlock 2:anxiety_c             -1.483e-02  1.827e-02  1.198e+02  -0.812
## blockBlock 1:stai_trait            -3.621e-04  4.344e-04  1.041e+03  -0.834
## blockBlock 2:stai_trait             1.199e-04  4.841e-04  1.189e+02   0.248
##                                    Pr(>|t|)    
## (Intercept)                          <2e-16 ***
## c_age_y_centered                     0.0550 .  
## c_sexF                               0.7790    
## cbq_fear                             0.2305    
## blockBlock 1                         0.0904 .  
## blockBlock 2                         0.0310 *  
## cbq_effortful_control                0.0422 *  
## anxiety_c                            0.7198    
## stai_trait                           0.2138    
## cbq_fear:blockBlock 1                0.3336    
## cbq_fear:blockBlock 2                0.0116 *  
## blockBlock 1:cbq_effortful_control   0.0457 *  
## blockBlock 2:cbq_effortful_control   0.0354 *  
## blockBlock 1:anxiety_c               0.6076    
## blockBlock 2:anxiety_c               0.4185    
## blockBlock 1:stai_trait              0.4047    
## blockBlock 2:stai_trait              0.8048    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## 
## Correlation matrix not shown by default, as p = 17 > 12.
## Use print(x, correlation=TRUE)  or
##     vcov(x)        if you need it

## optimizer (bobyqa) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')

# emm = emmeans(model3_4, ~ cbq_effortful_control*block )
# pairs(emm)
# # or for simple comparisons
# #pairs(emm, simple = "each")

anova(model3_4b, type = "III")

report(model3_4b)

## boundary (singular) fit: see help('isSingular')

## Random effect variances not available. Returned R2 does not account for random effects.

## boundary (singular) fit: see help('isSingular')

## Random effect variances not available. Returned R2 does not account for random effects.

## We fitted a linear mixed model (estimated using REML and BOBYQA optimizer) to
## predict DLPFC_avg_score with c_age_y_centered, c_sex, cbq_fear, block,
## cbq_effortful_control, anxiety_c and stai_trait (formula: DLPFC_avg_score ~
## c_age_y_centered + c_sex + cbq_fear * block + cbq_effortful_control * block +
## anxiety_c * block + stai_trait * block). The model included block as random
## effects (formula: list(~1 + block | record_id, ~1 | trial)). The model's
## explanatory power related to the fixed effects alone (marginal R2) is 0.06. The
## model's intercept, corresponding to c_age_y_centered = 0, c_sex = M, cbq_fear =
## 0, block = Block 3, cbq_effortful_control = 0, anxiety_c = 0 and stai_trait =
## 0, is at 0.38 (95% CI [0.33, 0.43], t(1229) = 14.69, p < .001). Within this
## model:
## 
##   - The effect of c age y centered is statistically non-significant and positive
## (beta = 2.66e-03, 95% CI [-4.82e-05, 5.36e-03], t(1229) = 1.93, p = 0.054; Std.
## beta = 0.06, 95% CI [-1.05e-03, 0.12])
##   - The effect of c sex [F] is statistically non-significant and negative (beta =
## -9.47e-04, 95% CI [-7.56e-03, 5.66e-03], t(1229) = -0.28, p = 0.779; Std. beta
## = -0.02, 95% CI [-0.14, 0.10])
##   - The effect of cbq fear is statistically non-significant and negative (beta =
## -3.37e-03, 95% CI [-8.89e-03, 2.14e-03], t(1229) = -1.20, p = 0.230; Std. beta
## = -0.06, 95% CI [-0.17, 0.04])
##   - The effect of block [Block 1] is statistically non-significant and negative
## (beta = -0.06, 95% CI [-0.13, 9.79e-03], t(1229) = -1.69, p = 0.090; Std. beta
## = 0.06, 95% CI [-0.07, 0.19])
##   - The effect of block [Block 2] is statistically significant and negative (beta
## = -0.09, 95% CI [-0.17, -8.98e-03], t(1229) = -2.18, p = 0.029; Std. beta =
## 0.44, 95% CI [0.30, 0.59])
##   - The effect of cbq effortful control is statistically significant and negative
## (beta = -8.50e-03, 95% CI [-0.02, -3.00e-04], t(1229) = -2.03, p = 0.042; Std.
## beta = -0.10, 95% CI [-0.19, -3.44e-03])
##   - The effect of anxiety c is statistically non-significant and positive (beta =
## 4.17e-03, 95% CI [-0.02, 0.03], t(1229) = 0.36, p = 0.720; Std. beta = 0.02,
## 95% CI [-0.09, 0.13])
##   - The effect of stai trait is statistically non-significant and positive (beta
## = 3.83e-04, 95% CI [-2.21e-04, 9.88e-04], t(1229) = 1.24, p = 0.214; Std. beta
## = 0.06, 95% CI [-0.04, 0.16])
##   - The effect of cbq fear × block [Block 1] is statistically non-significant and
## positive (beta = 3.81e-03, 95% CI [-3.91e-03, 0.01], t(1229) = 0.97, p = 0.334;
## Std. beta = 0.07, 95% CI [-0.07, 0.22])
##   - The effect of cbq fear × block [Block 2] is statistically significant and
## positive (beta = 0.01, 95% CI [2.63e-03, 0.02], t(1229) = 2.56, p = 0.011; Std.
## beta = 0.21, 95% CI [0.05, 0.38])
##   - The effect of block [Block 1] × cbq effortful control is statistically
## significant and positive (beta = 0.01, 95% CI [2.25e-04, 0.02], t(1229) = 2.00,
## p = 0.046; Std. beta = 0.13, 95% CI [2.58e-03, 0.26])
##   - The effect of block [Block 2] × cbq effortful control is statistically
## significant and positive (beta = 0.01, 95% CI [1.07e-03, 0.03], t(1229) = 2.13,
## p = 0.034; Std. beta = 0.16, 95% CI [0.01, 0.30])
##   - The effect of block [Block 1] × anxiety c is statistically non-significant
## and positive (beta = 8.46e-03, 95% CI [-0.02, 0.04], t(1229) = 0.51, p = 0.608;
## Std. beta = 0.04, 95% CI [-0.11, 0.19])
##   - The effect of block [Block 2] × anxiety c is statistically non-significant
## and negative (beta = -0.01, 95% CI [-0.05, 0.02], t(1229) = -0.81, p = 0.417;
## Std. beta = -0.07, 95% CI [-0.24, 0.10])
##   - The effect of block [Block 1] × stai trait is statistically non-significant
## and negative (beta = -3.62e-04, 95% CI [-1.21e-03, 4.90e-04], t(1229) = -0.83,
## p = 0.405; Std. beta = -0.06, 95% CI [-0.20, 0.08])
##   - The effect of block [Block 2] × stai trait is statistically non-significant
## and positive (beta = 1.20e-04, 95% CI [-8.30e-04, 1.07e-03], t(1229) = 0.25, p
## = 0.804; Std. beta = 0.02, 95% CI [-0.13, 0.17])
## 
## Standardized parameters were obtained by fitting the model on a standardized
## version of the dataset. 95% Confidence Intervals (CIs) and p-values were
## computed using a Wald t-distribution approximation.

Plot for interaction: Fear

interaction_eff3_4 <- ggpredict(model3_4, c("cbq_fear", "block"))
interaction_eff3_4

ggplot(interaction_eff3_4, aes(x, predicted, color = group)) +
  geom_line(linewidth = 1.2) +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high, fill = group), alpha = 0.2, color = NA) +
  theme_classic() +
  scale_color_manual(
    name = "Block Type",
    values = natparks.pals("Saguaro", 3),
    labels = c(
      "Block 1" = "Baseline",
      "Block 2" = "Stressful",
      "Block 3" = "Recovery"
    )
  ) +
  scale_fill_manual(
    name = "Block Type",
    values = natparks.pals("Saguaro", 3),
    labels = c(
      "Block 1" = "Baseline",
      "Block 2" = "Stressful",
      "Block 3" = "Recovery"
    )
  ) +
  scale_x_continuous(
    name = "Fear (CBQ scores)"
  ) +
  theme(
    axis.text = element_text(size = 19),
    axis.title = element_text(size = 20),
    legend.text = element_text(size = 19),
    legend.title = element_text(face = "bold", size = 19),
    legend.title.align = 0.5
  ) +
  scale_y_continuous(name = "DLPFC Neural Synchrony (predicted)", limits = c(0.29, .42), breaks = seq(0.29, .42, by = 0.02)) # Setting Y-axis limits and ticks

ggsave("figures/aim_3_F-Block_trait_DLPFC.png", width = 10, height = 8)

Plot for interaction: EC

interaction_eff3_4_b <- ggpredict(model3_4, c("cbq_effortful_control", "block"))
interaction_eff3_4_b

ggplot(interaction_eff3_4_b, aes(x, predicted, color = group)) +
  geom_line(linewidth = 1.2) +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high, fill = group), alpha = 0.2, color = NA) +
  theme_classic() +
  scale_color_manual(
    name = "Block Type",
    values = natparks.pals("Saguaro", 3),
    labels = c(
      "Block 1" = "Baseline",
      "Block 2" = "Stressful",
      "Block 3" = "Recovery"
    )
  ) +
  scale_fill_manual(
    name = "Block Type",
    values = natparks.pals("Saguaro", 3),
    labels = c(
      "Block 1" = "Baseline",
      "Block 2" = "Stressful",
      "Block 3" = "Recovery"
    )
  ) +
  scale_x_continuous(
    name = "Effortful Control (CBQ scores)"
  ) +
  theme(
    axis.text = element_text(size = 19),
    axis.title = element_text(size = 20),
    legend.text = element_text(size = 19),
    legend.title = element_text(face = "bold", size = 19),
    legend.title.align = 0.5
  ) +
  scale_y_continuous(name = "DLPFC Neural Synchrony (predicted)", limits = c(0.29, .42), breaks = seq(0.29, .42, by = 0.02)) # Setting Y-axis limits and ticks

ggsave("figures/aim_3_EC-Block_trait_DLPFC.png", width = 10, height = 8)

VLPFC

State - Model -block 2 reference

nirs_per_trial_plus_qst$block <- factor(nirs_per_trial_plus_qst$block, levels = c("Block 2", "Block 1", "Block 3"))

model3_5 <- lmer(
  VLPFC_avg_score ~ c_age_y_centered +
    c_sex + cbq_fear * block + cbq_effortful_control * block + anxiety_c * block + stai_state * block +
    (1 + block | record_id) + (1 | trial),
  data = nirs_per_trial_plus_qst,
  control = lmerControl(
    optimizer = "bobyqa",
    optCtrl = list(maxfun = 5e5)
  )
)

## boundary (singular) fit: see help('isSingular')

summary(model3_5)

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: VLPFC_avg_score ~ c_age_y_centered + c_sex + cbq_fear * block +  
##     cbq_effortful_control * block + anxiety_c * block + stai_state *  
##     block + (1 + block | record_id) + (1 | trial)
##    Data: nirs_per_trial_plus_qst
## Control: lmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 5e+05))
## 
## REML criterion at convergence: -4609.4
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.0446 -0.6848 -0.0171  0.6586  4.0821 
## 
## Random effects:
##  Groups    Name         Variance  Std.Dev. Corr       
##  record_id (Intercept)  8.226e-05 0.009070            
##            blockBlock 1 1.008e-04 0.010038 -0.98      
##            blockBlock 3 3.184e-05 0.005643 -0.86  0.94
##  trial     (Intercept)  0.000e+00 0.000000            
##  Residual               1.451e-03 0.038098            
## Number of obs: 1300, groups:  record_id, 109; trial, 4
## 
## Fixed effects:
##                                      Estimate Std. Error         df t value
## (Intercept)                         3.497e-01  2.073e-02  1.142e+02  16.865
## c_age_y_centered                   -7.064e-04  9.865e-04  1.229e+02  -0.716
## c_sexF                              1.868e-04  2.418e-03  1.236e+02   0.077
## cbq_fear                           -4.521e-03  2.174e-03  1.117e+02  -2.080
## blockBlock 1                        2.189e-02  2.803e-02  1.185e+02   0.781
## blockBlock 3                       -1.164e-02  2.683e-02  1.889e+02  -0.434
## cbq_effortful_control               5.676e-03  3.319e-03  1.167e+02   1.710
## anxiety_c                           1.915e-03  8.920e-03  1.117e+02   0.215
## stai_state                          1.694e-05  2.354e-04  1.118e+02   0.072
## cbq_fear:blockBlock 1               2.778e-03  2.959e-03  1.181e+02   0.939
## cbq_fear:blockBlock 3               6.517e-03  2.848e-03  1.916e+02   2.288
## blockBlock 1:cbq_effortful_control -9.672e-03  4.454e-03  1.187e+02  -2.172
## blockBlock 3:cbq_effortful_control -6.800e-03  4.265e-03  1.894e+02  -1.595
## blockBlock 1:anxiety_c              7.930e-04  1.216e-02  1.184e+02   0.065
## blockBlock 3:anxiety_c             -1.147e-02  1.163e-02  1.886e+02  -0.986
## blockBlock 1:stai_state            -1.493e-04  3.203e-04  1.179e+02  -0.466
## blockBlock 3:stai_state             1.546e-04  3.078e-04  1.903e+02   0.502
##                                    Pr(>|t|)    
## (Intercept)                          <2e-16 ***
## c_age_y_centered                     0.4753    
## c_sexF                               0.9385    
## cbq_fear                             0.0398 *  
## blockBlock 1                         0.4364    
## blockBlock 3                         0.6648    
## cbq_effortful_control                0.0899 .  
## anxiety_c                            0.8304    
## stai_state                           0.9428    
## cbq_fear:blockBlock 1                0.3497    
## cbq_fear:blockBlock 3                0.0232 *  
## blockBlock 1:cbq_effortful_control   0.0319 *  
## blockBlock 3:cbq_effortful_control   0.1125    
## blockBlock 1:anxiety_c               0.9481    
## blockBlock 3:anxiety_c               0.3253    
## blockBlock 1:stai_state              0.6420    
## blockBlock 3:stai_state              0.6161    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## 
## Correlation matrix not shown by default, as p = 17 > 12.
## Use print(x, correlation=TRUE)  or
##     vcov(x)        if you need it

## optimizer (bobyqa) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')

# emm = emmeans(model3_5, ~ cbq_effortful_control*block )
# pairs(emm)
# # or for simple comparisons
# #pairs(emm, simple = "each")

anova(model3_5, type = "III")

report(model3_5)

## boundary (singular) fit: see help('isSingular')

## Random effect variances not available. Returned R2 does not account for random effects.

## boundary (singular) fit: see help('isSingular')

## Random effect variances not available. Returned R2 does not account for random effects.

## We fitted a linear mixed model (estimated using REML and BOBYQA optimizer) to
## predict VLPFC_avg_score with c_age_y_centered, c_sex, cbq_fear, block,
## cbq_effortful_control, anxiety_c and stai_state (formula: VLPFC_avg_score ~
## c_age_y_centered + c_sex + cbq_fear * block + cbq_effortful_control * block +
## anxiety_c * block + stai_state * block). The model included block as random
## effects (formula: list(~1 + block | record_id, ~1 | trial)). The model's
## explanatory power related to the fixed effects alone (marginal R2) is 0.08. The
## model's intercept, corresponding to c_age_y_centered = 0, c_sex = M, cbq_fear =
## 0, block = Block 2, cbq_effortful_control = 0, anxiety_c = 0 and stai_state =
## 0, is at 0.35 (95% CI [0.31, 0.39], t(1275) = 16.86, p < .001). Within this
## model:
## 
##   - The effect of c age y centered is statistically non-significant and negative
## (beta = -7.06e-04, 95% CI [-2.64e-03, 1.23e-03], t(1275) = -0.72, p = 0.474;
## Std. beta = -0.02, 95% CI [-0.08, 0.04])
##   - The effect of c sex [F] is statistically non-significant and positive (beta =
## 1.87e-04, 95% CI [-4.56e-03, 4.93e-03], t(1275) = 0.08, p = 0.938; Std. beta =
## 4.69e-03, 95% CI [-0.11, 0.12])
##   - The effect of cbq fear is statistically significant and negative (beta =
## -4.52e-03, 95% CI [-8.78e-03, -2.56e-04], t(1275) = -2.08, p = 0.038; Std. beta
## = -0.12, 95% CI [-0.23, -6.77e-03])
##   - The effect of block [Block 1] is statistically non-significant and positive
## (beta = 0.02, 95% CI [-0.03, 0.08], t(1275) = 0.78, p = 0.435; Std. beta =
## -0.54, 95% CI [-0.68, -0.41])
##   - The effect of block [Block 3] is statistically non-significant and negative
## (beta = -0.01, 95% CI [-0.06, 0.04], t(1275) = -0.43, p = 0.664; Std. beta =
## -0.53, 95% CI [-0.66, -0.39])
##   - The effect of cbq effortful control is statistically non-significant and
## positive (beta = 5.68e-03, 95% CI [-8.35e-04, 0.01], t(1275) = 1.71, p = 0.087;
## Std. beta = 0.09, 95% CI [-0.01, 0.19])
##   - The effect of anxiety c is statistically non-significant and positive (beta =
## 1.92e-03, 95% CI [-0.02, 0.02], t(1275) = 0.21, p = 0.830; Std. beta = 0.01,
## 95% CI [-0.10, 0.13])
##   - The effect of stai state is statistically non-significant and positive (beta
## = 1.69e-05, 95% CI [-4.45e-04, 4.79e-04], t(1275) = 0.07, p = 0.943; Std. beta
## = 3.89e-03, 95% CI [-0.10, 0.11])
##   - The effect of cbq fear × block [Block 1] is statistically non-significant and
## positive (beta = 2.78e-03, 95% CI [-3.03e-03, 8.58e-03], t(1275) = 0.94, p =
## 0.348; Std. beta = 0.07, 95% CI [-0.08, 0.23])
##   - The effect of cbq fear × block [Block 3] is statistically significant and
## positive (beta = 6.52e-03, 95% CI [9.29e-04, 0.01], t(1275) = 2.29, p = 0.022;
## Std. beta = 0.17, 95% CI [0.02, 0.32])
##   - The effect of block [Block 1] × cbq effortful control is statistically
## significant and negative (beta = -9.67e-03, 95% CI [-0.02, -9.35e-04], t(1275)
## = -2.17, p = 0.030; Std. beta = -0.15, 95% CI [-0.29, -0.01])
##   - The effect of block [Block 3] × cbq effortful control is statistically
## non-significant and negative (beta = -6.80e-03, 95% CI [-0.02, 1.57e-03],
## t(1275) = -1.59, p = 0.111; Std. beta = -0.11, 95% CI [-0.24, 0.02])
##   - The effect of block [Block 1] × anxiety c is statistically non-significant
## and positive (beta = 7.93e-04, 95% CI [-0.02, 0.02], t(1275) = 0.07, p = 0.948;
## Std. beta = 5.27e-03, 95% CI [-0.15, 0.16])
##   - The effect of block [Block 3] × anxiety c is statistically non-significant
## and negative (beta = -0.01, 95% CI [-0.03, 0.01], t(1275) = -0.99, p = 0.324;
## Std. beta = -0.08, 95% CI [-0.23, 0.08])
##   - The effect of block [Block 1] × stai state is statistically non-significant
## and negative (beta = -1.49e-04, 95% CI [-7.78e-04, 4.79e-04], t(1275) = -0.47,
## p = 0.641; Std. beta = -0.03, 95% CI [-0.18, 0.11])
##   - The effect of block [Block 3] × stai state is statistically non-significant
## and positive (beta = 1.55e-04, 95% CI [-4.49e-04, 7.58e-04], t(1275) = 0.50, p
## = 0.616; Std. beta = 0.04, 95% CI [-0.10, 0.17])
## 
## Standardized parameters were obtained by fitting the model on a standardized
## version of the dataset. 95% Confidence Intervals (CIs) and p-values were
## computed using a Wald t-distribution approximation.

model3_5_x <- lmer(
  VLPFC_avg_score ~ c_age_y_centered +
    c_sex + cbq_fear * block + cbq_effortful_control * block + anxiety_c + stai_state * block +
    (1 + block | record_id) + (1 | trial),
  data = nirs_per_trial_plus_qst,
  control = lmerControl(
    optimizer = "bobyqa",
    optCtrl = list(maxfun = 5e5)
  )
)

## boundary (singular) fit: see help('isSingular')

model3_5_y <- lmer(
  VLPFC_avg_score ~ c_age_y_centered +
    c_sex + cbq_fear * block + cbq_effortful_control * block + anxiety_c + stai_state +
    (1 + block | record_id) + (1 | trial),
  data = nirs_per_trial_plus_qst,
  control = lmerControl(
    optimizer = "bobyqa",
    optCtrl = list(maxfun = 5e5)
  )
)

## boundary (singular) fit: see help('isSingular')

## Warning: Model failed to converge with 1 negative eigenvalue: -1.2e+01

performance::compare_performance(model3_5, model3_5_x, model3_5_y)

## Random effect variances not available. Returned R2 does not account for random effects.
## Random effect variances not available. Returned R2 does not account for random effects.
## Random effect variances not available. Returned R2 does not account for random effects.

## best model is model3_5 <-

Plot for Interaction: Fear

interaction_eff3_5 <- ggpredict(model3_5, c("cbq_fear", "block"))
interaction_eff3_5

ggplot(interaction_eff3_5, aes(x, predicted, color = group)) +
  geom_line(linewidth = 1.2) +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high, fill = group), alpha = 0.2, color = NA) +
  theme_classic() +
  scale_color_manual(
    name = "Block Type",
    values = natparks.pals("Saguaro", 3),
    labels = c(
      "Block 1" = "Baseline",
      "Block 2" = "Stressful",
      "Block 3" = "Recovery"
    )
  ) +
  scale_fill_manual(
    name = "Block Type",
    values = natparks.pals("Saguaro", 3),
    labels = c(
      "Block 1" = "Baseline",
      "Block 2" = "Stressful",
      "Block 3" = "Recovery"
    )
  ) +
  scale_x_continuous(
    name = "Fear (CBQ scores)"
  ) +
  theme(
    axis.text = element_text(size = 19),
    axis.title = element_text(size = 20),
    legend.text = element_text(size = 19),
    legend.title = element_text(face = "bold", size = 19),
    legend.title.align = 0.5
  ) +
  scale_y_continuous(name = "VLPFC Neural Synchrony (predicted)", limits = c(0.29, .42), breaks = seq(0.29, .42, by = 0.02)) # Setting Y-axis limits and ticks

ggsave("figures/aim_3_F-Block_state_VLPFC.png", width = 10, height = 8)

Plot for Interaction: EC

interaction_eff3_5_b <- ggpredict(model3_5, c("cbq_effortful_control", "block"))
interaction_eff3_5_b

ggplot(interaction_eff3_5_b, aes(x, predicted, color = group)) +
  geom_line(linewidth = 1.2) +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high, fill = group), alpha = 0.2, color = NA) +
  theme_classic() +
  scale_color_manual(
    name = "Block Type",
    values = natparks.pals("Saguaro", 3),
    labels = c(
      "Block 1" = "Baseline",
      "Block 2" = "Stressful",
      "Block 3" = "Recovery"
    )
  ) +
  scale_fill_manual(
    name = "Block Type",
    values = natparks.pals("Saguaro", 3),
    labels = c(
      "Block 1" = "Baseline",
      "Block 2" = "Stressful",
      "Block 3" = "Recovery"
    )
  ) +
  scale_x_continuous(
    name = "Effortful Control (CBQ scores)"
  ) +
  theme(
    axis.text = element_text(size = 19),
    axis.title = element_text(size = 20),
    legend.text = element_text(size = 19),
    legend.title = element_text(face = "bold", size = 19),
    legend.title.align = 0.5
  ) +
  scale_y_continuous(name = "VLPFC Neural Synchrony (predicted)", limits = c(0.29, .42), breaks = seq(0.29, .42, by = 0.02)) # Setting Y-axis limits and ticks

ggsave("figures/aim_3_EC-Block_state_VLPFC.png", width = 10, height = 8)

Trait - Model -block 2 reference

nirs_per_trial_plus_qst$block <- factor(nirs_per_trial_plus_qst$block, levels = c("Block 2", "Block 1", "Block 3"))

model3_6 <- lmer(
  VLPFC_avg_score ~ c_age_y_centered +
    c_sex + cbq_fear * block + cbq_effortful_control * block + anxiety_c * block + stai_trait * block +
    (1 + block | record_id) + (1 | trial),
  data = nirs_per_trial_plus_qst,
  control = lmerControl(
    optimizer = "bobyqa",
    optCtrl = list(maxfun = 5e5)
  )
)

## boundary (singular) fit: see help('isSingular')

summary(model3_6)

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: VLPFC_avg_score ~ c_age_y_centered + c_sex + cbq_fear * block +  
##     cbq_effortful_control * block + anxiety_c * block + stai_trait *  
##     block + (1 + block | record_id) + (1 | trial)
##    Data: nirs_per_trial_plus_qst
## Control: lmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 5e+05))
## 
## REML criterion at convergence: -4611.8
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.0562 -0.6751 -0.0129  0.6510  4.0949 
## 
## Random effects:
##  Groups    Name         Variance  Std.Dev. Corr       
##  record_id (Intercept)  7.708e-05 0.008780            
##            blockBlock 1 9.408e-05 0.009700 -0.98      
##            blockBlock 3 3.278e-05 0.005726 -0.88  0.95
##  trial     (Intercept)  0.000e+00 0.000000            
##  Residual               1.451e-03 0.038098            
## Number of obs: 1300, groups:  record_id, 109; trial, 4
## 
## Fixed effects:
##                                      Estimate Std. Error         df t value
## (Intercept)                         3.415e-01  2.039e-02  1.132e+02  16.751
## c_age_y_centered                   -7.246e-04  9.800e-04  1.209e+02  -0.739
## c_sexF                              9.351e-05  2.395e-03  1.214e+02   0.039
## cbq_fear                           -4.370e-03  2.163e-03  1.111e+02  -2.020
## blockBlock 1                        2.912e-02  2.765e-02  1.199e+02   1.053
## blockBlock 3                       -9.008e-03  2.659e-02  1.877e+02  -0.339
## cbq_effortful_control               5.794e-03  3.245e-03  1.158e+02   1.785
## anxiety_c                          -7.772e-04  9.070e-03  1.110e+02  -0.086
## stai_trait                          2.369e-04  2.407e-04  1.109e+02   0.984
## cbq_fear:blockBlock 1               2.490e-03  2.950e-03  1.194e+02   0.844
## cbq_fear:blockBlock 3               6.647e-03  2.850e-03  1.898e+02   2.332
## blockBlock 1:cbq_effortful_control -9.513e-03  4.366e-03  1.201e+02  -2.179
## blockBlock 3:cbq_effortful_control -7.135e-03  4.199e-03  1.879e+02  -1.699
## blockBlock 1:anxiety_c              3.699e-03  1.239e-02  1.197e+02   0.299
## blockBlock 3:anxiety_c             -1.119e-02  1.191e-02  1.873e+02  -0.940
## blockBlock 1:stai_trait            -3.566e-04  3.284e-04  1.192e+02  -1.086
## blockBlock 3:stai_trait             9.634e-05  3.161e-04  1.874e+02   0.305
##                                    Pr(>|t|)    
## (Intercept)                          <2e-16 ***
## c_age_y_centered                     0.4611    
## c_sexF                               0.9689    
## cbq_fear                             0.0457 *  
## blockBlock 1                         0.2944    
## blockBlock 3                         0.7352    
## cbq_effortful_control                0.0768 .  
## anxiety_c                            0.9319    
## stai_trait                           0.3271    
## cbq_fear:blockBlock 1                0.4003    
## cbq_fear:blockBlock 3                0.0207 *  
## blockBlock 1:cbq_effortful_control   0.0313 *  
## blockBlock 3:cbq_effortful_control   0.0909 .  
## blockBlock 1:anxiety_c               0.7657    
## blockBlock 3:anxiety_c               0.3486    
## blockBlock 1:stai_trait              0.2797    
## blockBlock 3:stai_trait              0.7609    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## 
## Correlation matrix not shown by default, as p = 17 > 12.
## Use print(x, correlation=TRUE)  or
##     vcov(x)        if you need it

## optimizer (bobyqa) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')

# emm = emmeans(model3_6, ~ cbq_effortful_control*block )
# pairs(emm)
# # or for simple comparisons
# #pairs(emm, simple = "each")

anova(model3_6, type = "III")

report(model3_6)

## boundary (singular) fit: see help('isSingular')

## Random effect variances not available. Returned R2 does not account for random effects.

## boundary (singular) fit: see help('isSingular')

## Random effect variances not available. Returned R2 does not account for random effects.

## We fitted a linear mixed model (estimated using REML and BOBYQA optimizer) to
## predict VLPFC_avg_score with c_age_y_centered, c_sex, cbq_fear, block,
## cbq_effortful_control, anxiety_c and stai_trait (formula: VLPFC_avg_score ~
## c_age_y_centered + c_sex + cbq_fear * block + cbq_effortful_control * block +
## anxiety_c * block + stai_trait * block). The model included block as random
## effects (formula: list(~1 + block | record_id, ~1 | trial)). The model's
## explanatory power related to the fixed effects alone (marginal R2) is 0.08. The
## model's intercept, corresponding to c_age_y_centered = 0, c_sex = M, cbq_fear =
## 0, block = Block 2, cbq_effortful_control = 0, anxiety_c = 0 and stai_trait =
## 0, is at 0.34 (95% CI [0.30, 0.38], t(1275) = 16.75, p < .001). Within this
## model:
## 
##   - The effect of c age y centered is statistically non-significant and negative
## (beta = -7.25e-04, 95% CI [-2.65e-03, 1.20e-03], t(1275) = -0.74, p = 0.460;
## Std. beta = -0.02, 95% CI [-0.08, 0.04])
##   - The effect of c sex [F] is statistically non-significant and positive (beta =
## 9.35e-05, 95% CI [-4.60e-03, 4.79e-03], t(1275) = 0.04, p = 0.969; Std. beta =
## 2.35e-03, 95% CI [-0.12, 0.12])
##   - The effect of cbq fear is statistically significant and negative (beta =
## -4.37e-03, 95% CI [-8.61e-03, -1.27e-04], t(1275) = -2.02, p = 0.044; Std. beta
## = -0.12, 95% CI [-0.23, -3.35e-03])
##   - The effect of block [Block 1] is statistically non-significant and positive
## (beta = 0.03, 95% CI [-0.03, 0.08], t(1275) = 1.05, p = 0.293; Std. beta =
## -0.54, 95% CI [-0.68, -0.41])
##   - The effect of block [Block 3] is statistically non-significant and negative
## (beta = -9.01e-03, 95% CI [-0.06, 0.04], t(1275) = -0.34, p = 0.735; Std. beta
## = -0.52, 95% CI [-0.66, -0.39])
##   - The effect of cbq effortful control is statistically non-significant and
## positive (beta = 5.79e-03, 95% CI [-5.72e-04, 0.01], t(1275) = 1.79, p = 0.074;
## Std. beta = 0.09, 95% CI [-9.07e-03, 0.19])
##   - The effect of anxiety c is statistically non-significant and negative (beta =
## -7.77e-04, 95% CI [-0.02, 0.02], t(1275) = -0.09, p = 0.932; Std. beta =
## -5.17e-03, 95% CI [-0.12, 0.11])
##   - The effect of stai trait is statistically non-significant and positive (beta
## = 2.37e-04, 95% CI [-2.35e-04, 7.09e-04], t(1275) = 0.98, p = 0.325; Std. beta
## = 0.05, 95% CI [-0.05, 0.16])
##   - The effect of cbq fear × block [Block 1] is statistically non-significant and
## positive (beta = 2.49e-03, 95% CI [-3.30e-03, 8.28e-03], t(1275) = 0.84, p =
## 0.399; Std. beta = 0.07, 95% CI [-0.09, 0.22])
##   - The effect of cbq fear × block [Block 3] is statistically significant and
## positive (beta = 6.65e-03, 95% CI [1.06e-03, 0.01], t(1275) = 2.33, p = 0.020;
## Std. beta = 0.18, 95% CI [0.03, 0.32])
##   - The effect of block [Block 1] × cbq effortful control is statistically
## significant and negative (beta = -9.51e-03, 95% CI [-0.02, -9.47e-04], t(1275)
## = -2.18, p = 0.030; Std. beta = -0.15, 95% CI [-0.29, -0.02])
##   - The effect of block [Block 3] × cbq effortful control is statistically
## non-significant and negative (beta = -7.14e-03, 95% CI [-0.02, 1.10e-03],
## t(1275) = -1.70, p = 0.089; Std. beta = -0.11, 95% CI [-0.24, 0.02])
##   - The effect of block [Block 1] × anxiety c is statistically non-significant
## and positive (beta = 3.70e-03, 95% CI [-0.02, 0.03], t(1275) = 0.30, p = 0.765;
## Std. beta = 0.02, 95% CI [-0.14, 0.19])
##   - The effect of block [Block 3] × anxiety c is statistically non-significant
## and negative (beta = -0.01, 95% CI [-0.03, 0.01], t(1275) = -0.94, p = 0.348;
## Std. beta = -0.07, 95% CI [-0.23, 0.08])
##   - The effect of block [Block 1] × stai trait is statistically non-significant
## and negative (beta = -3.57e-04, 95% CI [-1.00e-03, 2.88e-04], t(1275) = -1.09,
## p = 0.278; Std. beta = -0.08, 95% CI [-0.22, 0.06])
##   - The effect of block [Block 3] × stai trait is statistically non-significant
## and positive (beta = 9.63e-05, 95% CI [-5.24e-04, 7.17e-04], t(1275) = 0.30, p
## = 0.761; Std. beta = 0.02, 95% CI [-0.12, 0.16])
## 
## Standardized parameters were obtained by fitting the model on a standardized
## version of the dataset. 95% Confidence Intervals (CIs) and p-values were
## computed using a Wald t-distribution approximation.

performance::check_outliers(model3_6)

## OK: No outliers detected.
## - Based on the following method and threshold: cook (0.907).
## - For variable: (Whole model)

model3_6_x <- lmer(
  VLPFC_avg_score ~ c_age_y_centered +
    c_sex + cbq_fear * block + cbq_effortful_control * block + anxiety_c + stai_trait * block +
    (1 + block | record_id) + (1 | trial),
  data = nirs_per_trial_plus_qst,
  control = lmerControl(
    optimizer = "bobyqa",
    optCtrl = list(maxfun = 5e5)
  )
)

## boundary (singular) fit: see help('isSingular')

model3_6_y <- lmer(
  VLPFC_avg_score ~ c_age_y_centered +
    c_sex + cbq_fear * block + cbq_effortful_control * block + anxiety_c + stai_trait +
    (1 + block | record_id) + (1 | trial),
  data = nirs_per_trial_plus_qst,
  control = lmerControl(
    optimizer = "bobyqa",
    optCtrl = list(maxfun = 5e5)
  )
)

## boundary (singular) fit: see help('isSingular')

performance::compare_performance(model3_6, model3_6_x, model3_6_y)

## Random effect variances not available. Returned R2 does not account for random effects.
## Random effect variances not available. Returned R2 does not account for random effects.
## Random effect variances not available. Returned R2 does not account for random effects.

## the best model is model3_6 <-

Plot for Interaction: Fear

interaction_eff3_6 <- ggpredict(model3_6, c("cbq_fear", "block"))
interaction_eff3_6

ggplot(interaction_eff3_6, aes(x, predicted, color = group)) +
  geom_line(linewidth = 1.2) +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high, fill = group), alpha = 0.2, color = NA) +
  theme_classic() +
  scale_color_manual(
    name = "Block Type",
    values = natparks.pals("Saguaro", 3),
    labels = c(
      "Block 1" = "Baseline",
      "Block 2" = "Stressful",
      "Block 3" = "Recovery"
    )
  ) +
  scale_fill_manual(
    name = "Block Type",
    values = natparks.pals("Saguaro", 3),
    labels = c(
      "Block 1" = "Baseline",
      "Block 2" = "Stressful",
      "Block 3" = "Recovery"
    )
  ) +
  scale_x_continuous(
    name = "Fear (CBQ scores)"
  ) +
  theme(
    axis.text = element_text(size = 19),
    axis.title = element_text(size = 20),
    legend.text = element_text(size = 19),
    legend.title = element_text(face = "bold", size = 19),
    legend.title.align = 0.5
  ) +
  scale_y_continuous(name = "VLPFC Neural Synchrony (predicted)", limits = c(0.29, .42), breaks = seq(0.29, .42, by = 0.02)) # Setting Y-axis limits and ticks

ggsave("figures/aim_3_F-Block_trait_VLPFC.png", width = 10, height = 8)

Plot for Interaction: EC

interaction_eff3_6_b <- ggpredict(model3_6, c("cbq_effortful_control", "block"))
interaction_eff3_6_b

ggplot(interaction_eff3_6_b, aes(x, predicted, color = group)) +
  geom_line(linewidth = 1.2) +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high, fill = group), alpha = 0.2, color = NA) +
  theme_classic() +
  scale_color_manual(
    name = "Block Type",
    values = natparks.pals("Saguaro", 3),
    labels = c(
      "Block 1" = "Baseline",
      "Block 2" = "Stressful",
      "Block 3" = "Recovery"
    )
  ) +
  scale_fill_manual(
    name = "Block Type",
    values = natparks.pals("Saguaro", 3),
    labels = c(
      "Block 1" = "Baseline",
      "Block 2" = "Stressful",
      "Block 3" = "Recovery"
    )
  ) +
  scale_x_continuous(
    name = "Effortful Control (CBQ scores)"
  ) +
  theme(
    axis.text = element_text(size = 19),
    axis.title = element_text(size = 20),
    legend.text = element_text(size = 19),
    legend.title = element_text(face = "bold", size = 19),
    legend.title.align = 0.5
  ) +
  scale_y_continuous(name = "VLPFC Neural Synchrony (predicted)", limits = c(0.29, .42), breaks = seq(0.29, .42, by = 0.02)) # Setting Y-axis limits and ticks

ggsave("figures/aim_3_EC-Block_trait_VLPFC.png", width = 10, height = 8)

Graph: Anxiety Children vs. State Anxiety Parent

nirs_per_trial_plus_qst %>%
  distinct(record_id, .keep_all = TRUE) %>%
  select(record_id, DLPFC_avg_score, anxiety_c, stai_state, stai_trait) %>%
  ggplot(., aes(anxiety_c, DLPFC_avg_score)) +
  geom_point(size = 3) +
  labs(x = "Caregiver State Anxiety", y = "Dyadic Neural Synchrony Values") +
  geom_smooth(method = "lm") +
  scale_fill_brewer(palette = "Dark2") +
  scale_colour_brewer(palette = "Dark2") +
  guides(fill = FALSE, col = FALSE) +
  theme_apa() +
  theme(text = element_text(size = 25))

## `geom_smooth()` using formula = 'y ~ x'

## Warning: Removed 4 rows containing non-finite outside the scale range
## (`stat_smooth()`).

## Warning: Removed 4 rows containing missing values or values outside the scale range
## (`geom_point()`).

ggsave("figures/aim_3_corr_child_parent_state_anx.png", width = 12, height = 8)

## `geom_smooth()` using formula = 'y ~ x'

## Warning: Removed 4 rows containing non-finite outside the scale range (`stat_smooth()`).
## Removed 4 rows containing missing values or values outside the scale range
## (`geom_point()`).

Graph: Anxiety Children vs. State Anxiety Parent

nirs_per_trial_plus_qst %>%
  distinct(record_id, .keep_all = TRUE) %>%
  select(record_id, DLPFC_avg_score, anxiety_c, stai_state, stai_trait) %>%
  ggplot(., aes(stai_state, anxiety_c)) +
  geom_point(size = 3) +
  labs(x = "Caregiver State Anxiety", y = "Child Anxiety") +
  geom_smooth(method = "lm") +
  scale_fill_brewer(palette = "Dark2") +
  scale_colour_brewer(palette = "Dark2") +
  guides(fill = FALSE, col = FALSE) +
  theme_apa() +
  theme(text = element_text(size = 25))

## `geom_smooth()` using formula = 'y ~ x'

ggsave("figures/aim_3_corr_child_parent_state_anx.png", width = 12, height = 8)

## `geom_smooth()` using formula = 'y ~ x'

Graph: Anxiety Children vs. Trait Anxiety Parent

nirs_per_trial_plus_qst %>%
  distinct(record_id, .keep_all = TRUE) %>%
  select(record_id, DLPFC_avg_score, anxiety_c, stai_state, stai_trait) %>%
  ggplot(., aes(stai_trait, anxiety_c)) +
  geom_point(size = 3) +
  labs(x = "Caregiver Trait Anxiety", y = "Child Anxiety Score") +
  geom_smooth(method = "lm") +
  scale_fill_brewer(palette = "Dark2") +
  scale_colour_brewer(palette = "Dark2") +
  guides(fill = FALSE, col = FALSE) +
  theme_apa() +
  theme(text = element_text(size = 25))

## `geom_smooth()` using formula = 'y ~ x'

ggsave("figures/aim_3_corr_child_parent_trait_anx.png", width = 12, height = 8)

## `geom_smooth()` using formula = 'y ~ x'

Report-State

tab_model(model3_1, model3_3, model3_5, p.val = "kr", show.df = TRUE, collapse.ci = TRUE, p.style = "stars", digits = 3)

	avg score		DLPFC avg score		VLPFC avg score
Predictors	Estimates	df	Estimates	df	Estimates	df
(Intercept)	0.332 *** (0.290 – 0.373)	104.963	0.306 *** (0.242 – 0.370)	101.355	0.350 *** (0.309 – 0.391)	104.458
c age y centered	0.000 (-0.001 – 0.002)	101.804	0.003 * (0.000 – 0.005)	102.974	-0.001 (-0.003 – 0.001)	101.680
c sex [F]	0.000 (-0.004 – 0.004)	101.912	-0.001 (-0.008 – 0.006)	103.346	0.000 (-0.005 – 0.005)	101.853
cbq fear	-0.000 (-0.004 – 0.004)	103.468	0.008 * (0.001 – 0.014)	101.128	-0.005 * (-0.009 – -0.000)	102.240
block [Block 1]	0.024 (-0.033 – 0.080)	103.845	0.015 (-0.068 – 0.098)	99.690	0.022 (-0.034 – 0.077)	103.496
block [Block 3]	0.024 (-0.026 – 0.075)	103.539	0.077 (-0.009 – 0.163)	101.252	-0.012 (-0.065 – 0.042)	103.130
cbq effortful control	0.006 (-0.001 – 0.013)	106.772	0.005 (-0.005 – 0.015)	103.152	0.006 (-0.001 – 0.012)	106.685
anxiety c	-0.002 (-0.020 – 0.016)	103.256	-0.007 (-0.035 – 0.021)	101.692	0.002 (-0.016 – 0.020)	102.220
stai state	0.000 (-0.000 – 0.001)	103.500	0.000 (-0.001 – 0.001)	100.288	0.000 (-0.000 – 0.000)	102.285
cbq fear × block [Block 1]	-0.001 (-0.007 – 0.005)	103.783	-0.007 (-0.016 – 0.002)	100.614	0.003 (-0.003 – 0.009)	103.103
cbq fear × block [Block 3]	0.000 (-0.005 – 0.005)	104.840	-0.011 * (-0.021 – -0.002)	104.354	0.007 * (0.001 – 0.012)	104.223
block [Block 1] × cbq effortful control	-0.008 (-0.017 – 0.001)	103.963	-0.002 (-0.015 – 0.011)	100.069	-0.010 * (-0.019 – -0.001)	103.670
block [Block 3] × cbq effortful control	-0.009 * (-0.017 – -0.001)	103.726	-0.013 (-0.027 – 0.001)	102.091	-0.007 (-0.015 – 0.002)	103.359
block [Block 1] × anxiety c	0.005 (-0.020 – 0.029)	103.910	0.019 (-0.017 – 0.056)	102.225	0.001 (-0.023 – 0.025)	103.395
block [Block 3] × anxiety c	-0.003 (-0.025 – 0.019)	103.536	0.012 (-0.025 – 0.050)	102.577	-0.011 (-0.035 – 0.012)	102.934
block [Block 1] × stai state	-0.000 (-0.001 – 0.001)	103.698	-0.000 (-0.001 – 0.001)	99.742	-0.000 (-0.001 – 0.000)	102.986
block [Block 3] × stai state	0.000 (-0.000 – 0.001)	104.315	0.000 (-0.001 – 0.001)	101.561	0.000 (-0.000 – 0.001)	103.640
Random Effects
σ2	0.00		0.00		0.00
τ00	0.00 record_id		0.00 record_id		0.00 record_id
	0.00 trial		0.00 trial		0.00 trial
τ11	0.00 record_id.blockBlock 1		0.00 record_id.blockBlock 1		0.00 record_id.blockBlock 1
	0.00 record_id.blockBlock 3		0.00 record_id.blockBlock 3		0.00 record_id.blockBlock 3
ρ01	-1.00 record_id.blockBlock 1		-1.00 record_id.blockBlock 1		-0.98 record_id.blockBlock 1
	-1.00 record_id.blockBlock 3		-1.00 record_id.blockBlock 3		-0.86 record_id.blockBlock 3
N	109 record_id		109 record_id		109 record_id
	4 trial		4 trial		4 trial
Observations	1303		1254		1300
Marginal R2 / Conditional R2	0.090 / NA		0.058 / NA		0.076 / NA
p<0.05   ** p<0.01   *** p<0.001

Report-Trait

tab_model(model3_2, model3_4, model3_6, p.val = "kr", show.df = TRUE, collapse.ci = TRUE, p.style = "stars", digits = 3)

	avg score		DLPFC avg score		VLPFC avg score
Predictors	Estimates	df	Estimates	df	Estimates	df
(Intercept)	0.323 *** (0.283 – 0.364)	104.531	0.295 *** (0.232 – 0.357)	101.186	0.342 *** (0.301 – 0.382)	104.122
c age y centered	0.000 (-0.001 – 0.002)	101.819	0.003 (-0.000 – 0.005)	102.855	-0.001 (-0.003 – 0.001)	101.686
c sex [F]	0.000 (-0.004 – 0.004)	101.838	-0.001 (-0.007 – 0.006)	103.108	0.000 (-0.005 – 0.005)	101.787
cbq fear	0.000 (-0.004 – 0.004)	103.479	0.008 * (0.001 – 0.015)	101.013	-0.004 * (-0.009 – -0.000)	102.251
block [Block 1]	0.033 (-0.023 – 0.088)	103.818	0.027 (-0.055 – 0.108)	99.957	0.029 (-0.026 – 0.084)	103.588
block [Block 3]	0.031 (-0.019 – 0.081)	103.482	0.089 * (0.004 – 0.174)	101.645	-0.009 (-0.062 – 0.044)	103.208
cbq effortful control	0.006 (-0.001 – 0.012)	106.529	0.005 (-0.005 – 0.015)	102.954	0.006 (-0.001 – 0.012)	106.449
anxiety c	-0.005 (-0.023 – 0.013)	103.217	-0.011 (-0.039 – 0.017)	101.303	-0.001 (-0.019 – 0.017)	102.137
stai trait	0.000 (-0.000 – 0.001)	103.314	0.001 (-0.000 – 0.001)	100.349	0.000 (-0.000 – 0.001)	102.088
cbq fear × block [Block 1]	-0.001 (-0.007 – 0.005)	103.779	-0.007 (-0.016 – 0.001)	100.347	0.002 (-0.003 – 0.008)	103.095
cbq fear × block [Block 3]	0.000 (-0.005 – 0.005)	104.645	-0.011 * (-0.021 – -0.002)	104.197	0.007 * (0.001 – 0.012)	104.013
block [Block 1] × cbq effortful control	-0.008 (-0.017 – 0.001)	103.962	-0.002 (-0.015 – 0.011)	100.246	-0.010 * (-0.018 – -0.001)	103.711
block [Block 3] × cbq effortful control	-0.010 * (-0.018 – -0.002)	103.605	-0.014 * (-0.027 – -0.000)	102.089	-0.007 (-0.015 – 0.001)	103.277
block [Block 1] × anxiety c	0.008 (-0.017 – 0.033)	103.911	0.024 (-0.013 – 0.060)	101.559	0.004 (-0.021 – 0.028)	103.373
block [Block 3] × anxiety c	-0.002 (-0.024 – 0.021)	103.572	0.015 (-0.023 – 0.054)	102.556	-0.011 (-0.035 – 0.012)	102.956
block [Block 1] × stai trait	-0.000 (-0.001 – 0.000)	103.689	-0.000 (-0.001 – 0.000)	99.666	-0.000 (-0.001 – 0.000)	102.986
block [Block 3] × stai trait	-0.000 (-0.001 – 0.001)	103.685	-0.000 (-0.001 – 0.001)	101.881	0.000 (-0.001 – 0.001)	102.967
Random Effects
σ2	0.00		0.00		0.00
τ00	0.00 record_id		0.00 record_id		0.00 record_id
	0.00 trial		0.00 trial		0.00 trial
τ11	0.00 record_id.blockBlock 1		0.00 record_id.blockBlock 1		0.00 record_id.blockBlock 1
	0.00 record_id.blockBlock 3		0.00 record_id.blockBlock 3		0.00 record_id.blockBlock 3
ρ01	-1.00 record_id.blockBlock 1		-1.00 record_id.blockBlock 1		-0.98 record_id.blockBlock 1
	-1.00 record_id.blockBlock 3		-1.00 record_id.blockBlock 3		-0.88 record_id.blockBlock 3
N	109 record_id		109 record_id		109 record_id
	4 trial		4 trial		4 trial
Observations	1303		1254		1300
Marginal R2 / Conditional R2	0.092 / NA		0.060 / NA		0.077 / NA
p<0.05   ** p<0.01   *** p<0.001

4. Correlations with all

report(cor.test(block1$stai_state, block1$cbq_fear, method = "pearson", conf.level = 0.95))

## Effect sizes were labelled following Funder's (2019) recommendations.
## 
## The Pearson's product-moment correlation between block1$stai_state and
## block1$cbq_fear is positive, statistically significant, and small (r = 0.18,
## 95% CI [0.08, 0.27], t(432) = 3.74, p < .001)

report(cor.test(block1$stai_trait, block1$cbq_fear, method = "pearson", conf.level = 0.95))

## Effect sizes were labelled following Funder's (2019) recommendations.
## 
## The Pearson's product-moment correlation between block1$stai_trait and
## block1$cbq_fear is positive, statistically significant, and very small (r =
## 0.09, 95% CI [1.01e-04, 0.19], t(432) = 1.97, p = 0.050)

report(cor.test(block1$anxiety_c, block1$cbq_fear, method = "pearson", conf.level = 0.95))

## Effect sizes were labelled following Funder's (2019) recommendations.
## 
## The Pearson's product-moment correlation between block1$anxiety_c and
## block1$cbq_fear is positive, statistically significant, and very large (r =
## 0.46, 95% CI [0.39, 0.54], t(432) = 10.90, p < .001)

report(cor.test(block1$stai_state, block1$cbq_effortful_control, method = "pearson", conf.level = 0.95))

## Effect sizes were labelled following Funder's (2019) recommendations.
## 
## The Pearson's product-moment correlation between block1$stai_state and
## block1$cbq_effortful_control is negative, statistically significant, and small
## (r = -0.16, 95% CI [-0.25, -0.07], t(432) = -3.43, p < .001)

report(cor.test(block1$stai_trait, block1$cbq_effortful_control, method = "pearson", conf.level = 0.95))

## Effect sizes were labelled following Funder's (2019) recommendations.
## 
## The Pearson's product-moment correlation between block1$stai_trait and
## block1$cbq_effortful_control is negative, statistically not significant, and
## tiny (r = -0.03, 95% CI [-0.12, 0.06], t(432) = -0.62, p = 0.536)

report(cor.test(block1$anxiety_c, block1$cbq_effortful_control, method = "pearson", conf.level = 0.95))

## Effect sizes were labelled following Funder's (2019) recommendations.
## 
## The Pearson's product-moment correlation between block1$anxiety_c and
## block1$cbq_effortful_control is positive, statistically not significant, and
## tiny (r = 0.04, 95% CI [-0.05, 0.13], t(432) = 0.82, p = 0.413)

report(cor.test(block1$stai_state, block1$stai_trait, method = "pearson", conf.level = 0.95))

## Effect sizes were labelled following Funder's (2019) recommendations.
## 
## The Pearson's product-moment correlation between block1$stai_state and
## block1$stai_trait is positive, statistically significant, and very large (r =
## 0.78, 95% CI [0.75, 0.82], t(432) = 26.31, p < .001)

report(cor.test(block1$stai_trait, block1$anxiety_c, method = "pearson", conf.level = 0.95))

## Effect sizes were labelled following Funder's (2019) recommendations.
## 
## The Pearson's product-moment correlation between block1$stai_trait and
## block1$anxiety_c is positive, statistically significant, and large (r = 0.32,
## 95% CI [0.24, 0.41], t(432) = 7.10, p < .001)

report(cor.test(block1$stai_state, block1$anxiety_c, method = "pearson", conf.level = 0.95))

## Effect sizes were labelled following Funder's (2019) recommendations.
## 
## The Pearson's product-moment correlation between block1$stai_state and
## block1$anxiety_c is positive, statistically significant, and medium (r = 0.29,
## 95% CI [0.20, 0.37], t(432) = 6.22, p < .001)

report(cor.test(block1$cbq_effortful_control, block1$anxiety_c, method = "pearson", conf.level = 0.95))

## Effect sizes were labelled following Funder's (2019) recommendations.
## 
## The Pearson's product-moment correlation between block1$cbq_effortful_control
## and block1$anxiety_c is positive, statistically not significant, and tiny (r =
## 0.04, 95% CI [-0.05, 0.13], t(432) = 0.82, p = 0.413)